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ABSTRACT 

The  main  goal  of  thus  research  is  to  model  a  flexible  missile  with  structural 
flexibility  utilizing  the  Equivalent  Rigid  Link  System  (ERLS)  with  an  enhanced 
natural  mode  discretization.  Dynamic  analysis  of  the  flexible  missile  in 
2- Dimension  motion  is  presented. 

Computer  simulation  is  performed  where  the  pitch  angle  of  the  missile  is 
controlled  with  a  ngid*body  controller.  The  effects  of  increasing  payloads  and 
speed  to  the  system  performance  are  discussed. 
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THESIS  DISCLAIMER 


The  reader  is  cautioned  that  computer  programs  developed  in  this  research 
may  not  have  been  exercised  for  all  cases  of  interest.  While  every  effort  has  been 
made,  within  the  time  available,  to  ensure  that  the  programs  are  free  of  compu¬ 
tational  and  logic  errors,  they  cannot  be  considered  validated.  Any  application 
of  these  programs  without  additional  verification  is  at  the  risk  of  the  user. 
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I.  INTRODUCTION 


A.  BACKGROUND 

With  the  introduction  of  long  slender  missiles  such  as  the  Vanguard,  the 
Redstone,  and  various  ballistic  missiles,  the  problem  of  the  structural  flexibility 
became  severe  [Ref.  1].  Due  to  the  limited  thrust  available  from  rocket  engines, 
these  missiles  had  to  be  as  light  as  possible.This  meant  a  sacrifice  in  structural 
rigidity. 

Missile  flexure  causes  additional  aerodynamic  loads  which  in  turn  cause  ad¬ 
ditional  flexure.  Coupling  occurs  between  the  elastic  modes  and  the  control  sys¬ 
tem  as  the  control  system  gyros  sense  the  flexure  motion  and  the  rigid  body 
motion.  It  has  become  necessary  to  actively  control  the  flexible  structure  and 
thereby  reduce  the  structural  loads  and  improve  the  vehicle  response  such  as  po¬ 
sition,  velocity  and  acceleration.  Reduced  structural  loads  will  also  offer  potential 
for  reduced  bending  stress  and  fatigue  problems. 

The  objective  of  this  thesis  is  to  develop  a  dynamic  model  for  a  flexible  missile 
and  study  the  dynamic  behavior  of  the  flexible  missile.  Simulation  is  a  valuable 
tool  in  the  design  of  new  missile  systems  and  in  the  modification  or  evaluation 
of  existing  systems.  A  missile  simulation  allows  the  engineer  to  evaluate  his  design 
without  the  expense  of  actually  building  and  Hying  the  actual  missile.  System 
dynamics  can  be  investigated  through  simulation  with  a  substantial  savings  in 
time  and  expense  [Ref.  2]. 

B.  LITERATURE  REVIEW 

Flexible  missile  modeling  centers  on  the  relationship  between  the  large,  rigid- 
body  motion  and  the  small  motions  due  to  structural  flexibility. 

Jenkins  [Ref.  2]  expresses  some  techniques  used  in  deriving  the  equations  of 
motion  of  a  rigid  missile  for  a  six  degree-of-freedom  (6-DOF)  simulation.  The 
rigid  missiles  are  characterized  by  their  larger  size  and  low  ratio  of  payload  to 
total  weight.  Rigid  missile  dynamic  equations  were  developed  using  the 
Newton-Euler  Method.  The  moments  and  forces  along  with  the  mass  and  mo- 


merits  of  inertia  are  assumed  to  be  known  in  the  body  coordinate  frame.  The 
transformation  between  global  and  local  coordinate  frame  is  achieved  with  a 
non-homogenous  coordinate  transformation  matrix  [Ref.  3:  pp.  342].  The  result¬ 
ing  rigid-body  equations  of  motion  produced  the  understanding  for  the  derivation 
of  the  dynamic  equations  of  the  flexible  missile. 

The  Equivrdent  Rigid  Link  System  [Ref.  4]  describes  the  large-motion 
kinematics  of  the  system  by  the  ERLS  motion  and  the  small  motion  relative  to 
the  ERLS.  The  application  of  the  finite  element  techniques  and  Lagrangian  dy¬ 
namics  produces  two  sets  of  coupled,  nonlinear,  ordinary  differential  equations 
of  motion,  of  which  one  set  is  for  the  large  motions  and  the  other  set  for  the  small 
motions.  The  small  motion  is  described  by  the  superposition  of  vibration  modes. 
The  modes  of  the  vibration  of  the  flexible  bar  was  derived  with  the  simple-beam 
theory  [Ref.  5:  pp.221j.  In  simple  beam  theory,  it  is  assumed  that  the  rotation  of 
the  element  is  insignificant  compared  to  the  vertical  translation  and  the  shear 
deformation  is  small  relative  to  the  bending  deformation.  This  assumption  is  valid 
if  the  ratio  between  the  length  of  the  bar  and  its  height  is  relatively  large  (more 
than  10).  The  set  of  large  motion  equations  are  non-linear  in  both  the  large  and 
small  motion  variables  while  the  set  of  small  motion  equations  are  linear  in  the 
small  motion  variable  and  non-linear  in  the  large  motion  variables.  A  solution 
technique  called  the  Sequential  Integration  Method  [Ref.  6]  was  developed  which 
allows  efficient  simulations  of  systems  with  inertia  coupled  motions  hav  ing  non¬ 
linear  slow  motion  (large  motion)  with  linear  fast  motion  (small  motion).  The 
ERLS  model  presents  a  complete,  efficient  dynamic  model  able  to  describe  large 
motion,  small  motion  and  their  coupling. 

An  ERLS  model  of  a  flexible  spacecraft  boom  was  developed  and  a  computer 
simulation  was  performed  [Ref.  7].  The  equations  of  motion  were  solved  using 
the  Sequential  Integration  Method.  A  spatial  finite  discretization  of  the  boom 
structure  and  the  application  of  an  assumed  polynomial  modal  response  were 
utilize^  in  the  approximate  solution  to  the  equations  of  motion.  This  work  was 
the  basis  work  for  the  modeling  of  the  flexible  missiles. 


Ganon  [Ref.  8]  performed  an  experimental  validation  of  the  ERLS  dynamic 
mode!  in  a  vertical  plane  motion.  The  small  motion  was  modeled  using  a  shape 
matrix  derived  from  superposition  of  natural  modes.  The  agreement  between  the 
simulation  results  and  the  experimental  results  justify  the  application  of  the 
ERLS  using  a  natural-mode  shape  matrix  to  model  the  dynamic  response  of 
flexiole  missile. 

C.  OUTLINE 

In  this  study,  the  ERLS  dynamic  model  is  used  to  derive  the  system  equations 
of  motion  for  a  flexible  missile  in  2-D  motion.  Dynamic  response  is  predicted  by 
solving  the  equations  of  motion  using  the  Sequential  Integration  Method.  The 
application  of  an  assumed  natural  mode  shape  and  the  spatial  finite  element 
discretization  of  the  missile  provide  an  approximate  solution.  Computer  simu¬ 
lation  for  the  flexible  missile  is  performed  using  MATLAB  on  the  MACINTOSH 
computer.  A  rigid  body  controller  is  included  in  the  simulation  to  control  the 
pitch  angle  of  the  flexible  missile. 

The  ERLS  dynamic  model  of  the  flexible  missile  is  presented  in  Chapter  Two. 
A  rigid  body  controller  for  the  flexible  missile  is  described  in  Chapter  Three.  The 
computer  simulation  methodology  and  simulation  results  are  presented  in  Chap¬ 
ter  Four.  The  conclusions  and  recommendations  for  future  work  are  presented 
in  Chapter  Five. 
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II.  MODEL  FORMULATION  FOR  A  FLEXIBLE  MISSILE 


A,  KINEMATICS 

In  our  model,  2-D  inertial  reference  frame,  i.e.,  X  and  V,  is  used  for  the  global 
motion  as  shown  in  Fig.l.  The  body-fixed  coordinate  frame,  i.e.,  x  and  y,  is  se¬ 
lected  to  describe  the  missile  local  motion. 


Figure  1.  The  general  configuration  of  the  missile 


The  following  assumptions  were  made  for  the  flexible  missile: 

(1)  Material  density  is  constant  throughout  the  body  and  the  steel  was  chcscn 
for  modeling. 

(2)  A  uniform  circular  cross  section  is  assumed. 

The  geometric  configuration  parameters  and  material  properties  of  the  flexi¬ 
ble  missile  are  listed  as 
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Diameter  «■  0.1 2  meter 
Length  =*4.0  meter 


Material  Density  “7861.05  kg/mJ 
Young's  Modulus”2.0  x  10"  pascal 

The  concept  of  the  ERLS  is  applied  to  model  the  kinematics  of  flexible  mis¬ 
sile.  The  main  idea  of  the  Equivalent  Rigid  Link  System  (Fig.2)  is  to  separate  the 
kinematics  of  the  flexible  body  into  large  and  small  motions. 


Figure  2.  Equivalent  Rigid  Link  System  (ERLS) 
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The  large  motion  of  the  missile  will  be  represented  by  missile  base  position  in 
X  direction  XQy  missile  base  position  in  V  direction  T0,  and  missile  pitch  angle  6. 
The  small  motions  resulting  from  flexible  motion  are  measured  relative  to  the  lo¬ 
cal  coordinate  frame  x,  y.  v(Q)  and  <t>(0)  are  the  nodal  displacement  ana  slope 
of  the  missile  base  lespectively.  The  absolute  (global)  position  of  a  point  on  the 
flexible  missile  is  obtained  from  coordinate  transformations.  The  global  position 
(R)  of  any  point  position  can  be  defined  using  ERLS  in  terms  of  a  local  position 
vector  ( r ),  a  deformation  vector  (d)y  and  a  coordinate  transformation  matrix  (W), 
i.e., 

R^.lV{T  +  d)  (2-1) 


The  transformation  matrix  (W)  between  the  large  motion  and  small  motion  co¬ 
ordinate  is 


IV  = 


1 

Yq 


0 

cos($) 

sin(d) 


0 

sin(0) 

cos(d) 


(2  -  2) 


and  the  local  rigid  body  position  vector  is 


r  = 


(2-3) 


The  deformation  vector  (d)  is  expressed  in  terms  of  a  nodal  displacement 
vector  U  and  a  shape  function  <f>  as 


d  =  <pU  (2-4) 

where 


£?  =  [»( 0)  4>(0)]r 


(2  -  5) 
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The  shape  function  <f>  was  derived  utilizing  a  natural-mode  superposition  and 
a  finite  element  concept.  The  Finite  Element  Method  (FEM)  was  utilized  to 
discretize  the  flexible  body  displacements  and  assigning  the  nodes.  Displacements 
are  for  each  pcint  along  the  missile,  a  function  of  location  and  time,  ana  it  is 
necessary  to  discretize  the  deformations  to  obtain  an  ordinary  differential 
equation.  The  natural  mode  shape  function  of  a  beam  is  used  to  represent  the 
flexural  motion  of  the  flexible  missile.  Only  the  first  two  mode  shapes  are  used. 
The  flexible  missile  is  modeled  as  a  continuous  Euler-Bernoulli  free-free  beam, 
neglecting  shear  deformation  and  rotary  inertia  effects.  The  nodal  points  are  the 
base  points  of  the  flexible  missile.  Appendix  A  shows  the  mode  shape  function 
derivation.  In  this  case,  0  is  found  as, 


4>  = 


”0  0 
0  0 
a  b 


(2-6) 


where  a  and  b  are  defined  as  following, 
a  -  F^C^  cos  /fix  4-  cosh  p{x)  +  ( sin  p{x  +  sinh  p{x)) 

4'  F\(C2{  cos  p2x  +  c°sh  Pix)  +  (  sin  @2X  +  sinh  Pzx))  (2  “  7) 

b  —  cos  jJjX  4-  cosh  p{x)  +  ( sin  p{x  4  sinh  /lj..v)) 

4  F4(C2(  cos  p2x  4-  cosh  p2x)  4  ( sin  P2x  4-  sinh  P2x))  (2  -  8) 

Lagrangian  dynamics  is  used  to  acquire  the  equations  of  motion  and  the 
kinetic  energy  of  the  system  will  be  needed  in  the  development  of  the  Lagrangian 
expression.  The  absolute  velocity  is  listed  as  follows, 


R  =  W(f  4  d)  f  Wd 


(2  -  9) 
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B.  KINETICS 

Lagrangian  Dynamics  is  a  systematic  way  to  derive  equations  of  motion  for 
complex  systems  like  flexible  missiles.  Lagrange's  equation  is  written  as, 


d  ,  dKE  v  SKE  dPE 
dt  dqi  dq{  dqt 


(/  =  1,2 


(2  -  10) 


where 

K.E.  =  Kinetic  Energy 
P.E.  =  Potential  Energy 
<7,  =  Generalized  coordinates 
Q  =  Generalized  forces 
n  =  Number  of  degrees  of  freedom 
The  generalized  coordinates  are  chosen  to  be, 

q  =  [X0  Y„  e  v(0)  O(0)]r  (2-11) 

The  kinetic  energy  of  the  system  is  defined  as  follows, 


K.E.  =  -j\RTRdm  (2-12) 

By  substituting  Eq.(2-1)  into  Eq.  (2-12),  the  kinetic  energy  is  written  as 


KE  =  y{(f7‘lFrlf7  +  2rTWTW4>U  +  2rTWTW4>U 

+  UT4>TWTW4>U  +  2  UT<t>TWTW(f)U  +  UT4>TWr\V4>UT)dm  (2-13) 

The  potential  energy  of  the  system  includes  the  strain  energy  of  the  flexible 
missile  and  the  gravitational  potential  energy. 

P.E.  =  ±-\UTVTCVUdx -\RTgdm  (2-14) 
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where 


T  =  Spatial  derivative  of  the  shape  function 
C  -  Rigidity  Matrix 
g  =*  Gravitational  acceleration  vector 

Appendix  B  includes  the  development  of  the  equations  of  motions  using 
Lagrange  equations. 

Generalized  forces  will  be  found  by  virtual  work  principle.  It  is  assumed  that 
the  only  force  is  the  thrust  force  which  is  applied  to  the  base  of  the  flexible  missile 
as  shown  in  Fig. 3.  Other  forces  like  aerodynamic  and  damping  forces  are  neg¬ 
lected. 


Figure  3.  Applied  Forces  on  Flexible  Missile 


Tlfl^dirust  force  is  composed  of  two  components  i.e., 
Tx  =  rcos(<5  -F  <J>(0)) 


(2-15) 
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Ty~  7  sin(<5  +  0(0)) 


(2-16) 


and  c  moment  applied  on  the  missile  is, 

Me  =  -T cos(<5  4-  0(0))v(0)  (2-17) 

By  applying  virtual  work  principle,  the  generalized  forces  of  flexible  missile  is 
found  below\ 


7cos(0)  -  T5  sin (0)  -  7<D(0)  sin(0) 
7  sin(0)  +  TS  cos(0)  +70(0)  cos(0) 
-7v(0) 


(2-18) 


7<5  +  70(0) 
~7v(0) 


(2-19) 


where 

T  -  Force 

5  =  Control  Angle  being  assumed  small 
Fj  =  Large  motion  generalized  force  vector 
Fu  ~  Small  motion  generalized  force  vector 

Appendix  C  shows  the  derivation  of  the  generalized  forces  using  the  Virtual 
Work  Principle. 

The  derivation  of  the  equations  of  motion  from  the  Lagrangian  formulation 
yields  two  sets  of  coupled  equations.  One  set  describes  the  large  motions  and  the 
another  set  describes  the  small  motions.  These  two  sets  of  equations  are  non¬ 
linear,  coupled,  second-order,  ordinary  differential  equations  of  the  form, 


{  +  Mqn  U  =  Fq 


(2  -  20) 
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(2-21) 


Mnq  t  +  Mnn  U  +  GnU  +  KnU  =  /; 
where: 

Mqq  =  3x3  effective  inertia  matrix  for  large  motions 

Mqn  =  3x2  coupled  inertia  matrix  cf  the  small  motion  effect  on  large 

motions 

Fq  =  3x1  load  vector  for  the  large  motions 

M„q  s  2x3  coupled  inertia  matrix  of  the  large  motion  effect  on  small 
motions 

Mnn  =  2x2  effective  inertia  matrix  for  small  motions 

G„  s  2x2  gyroscopic  matrix 
K„  =  2x2  stiffness  matrix 

Fn  =  2x1  load  vector  for  the  small  motions 

£  =  3x1  vector,  generalized  coordinates  of  the  large  motions 

U  =  2x1  vector,  generalized  coordinates  of  the  small  motions 

The  detailed  development  of  the  equations  of  motion  and  definitions  of  the  terms 
are  described  in  Appendix  B. 
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HI.  THE  DEVELOPMENT  OF  A  RIGID-BODY  CONTROLLER 

In  our  research,  the  controller  Is  developed  based  on  the  rigid-body  assump¬ 
tion.  The  purpose  of  including  the  rigid-body  controller  is  to  perform  computer 
simulation  and  study  the  dynamic  behavior  of  the  flexible  missile.  The  pitch  an¬ 
gle  of  the  missile  is  controlled  by  the  rigid-body  controller  which  is  a  single  input 
single  output  (SISO)  controller.  A  general  configuration  of  the  flexible  missile 
and  rigid-body  controller  are  shown  in  Fig.4. 


(fy)  (£) 


Figure  4.  The  flexible  missile  »<*h  rigid-body  controller 


The  pitch  angle  is  the  control  variable  and  the  desired  states  are  desired  pitch 
angle  0d,  desired  angular  velocity  6d,  and  desired  angular  acceleration  6d.  The 
control  angle  (<5)  is  assumed  small. 

The  desired  error  function  is  defined  as 

£  +  Kvt  4-  Kpt  =  6  (3  -  1 ) 


where, 
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Kp  —  The  position  feedback  gain. 

Kv  ~  The  velocity  f  ,-edback  gain. 

£  =  6  -  dd  =  The  error  in  the  position. 

£  =  d  -  6d  =  The  error  in  velocity. 

£  =  8  -  Qd  —  The  error  in  acceleration. 

The  control  design  begins  with  the  equation  of  rigid-body  (large  motion)  that 
is  obtained  by  modifying  and  expanding  Eq.(2-20)  as 


— 1 

5 

j: 
_ 1 

'm' 

i 

i 

> 

0\ 

n 

+  +  O 

LA/2I 

i 

i 

lAl  M 

J>2. 

Where, 


Mqq(\A)  \fqq(  1,2) 

A*w(2,l)  ^w(2.2) 


*lqq\l>J) 

Mqq{  2,3) 


(3-3) 


(3-4) 


=  M„(3,2)]  (3-5) 

M22  -  [A/w(3,3)]  (3  —  6) 


^0 


(3  -  7) 


(3  -  X) 


f  =  Gravity,  centrifugal  and  coriohs  forces. 
f2  =  Gravity,  centrifugal  and  coriolis  forces. 
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-  \T  cos(0) 
h\  = 

[  T  sin(0) 

h2  =  [0] 

-  —  Tsin(0) 
1  =  T  cos(d) 


(3-9) 


(3  -  10) 


h2  =  [0] 


Eq.(3-2)  can  be  rewritten  in  a  tensor  form  as 
A/n^i  4-  A/12?/2  =./i  +  +  £q<5  (3  —  11) 

A^i^i  4-  M 22*1 2  ~fi  (3  —  12) 


From  Eq.(3-ll), 


rh  -  a/h1 

/,  4-  Aj  4-  />j<5  -  A/12t/2 

(3-13) 

Substituting  Eq.(3-13)  into  Eq.(3-12),  we  find, 

Arj 2  --F-B5  (3-14) 

Where, 

A  =  M22  - 

(3-15) 

f  =7: -W2] 

(3-16) 

B  *  A^A/f^ 

(3-17) 
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Eq.(3-1)  can  be  rewritten  as 


(9  -  ea  +  Kv(e  -  ed)  +  Kp(6  -  $d)  =  0  (3-18) 

By  substituting  Eq.(3-14)  into  Eq.(3-18),  we  find  the  control  angle  (b) 

S  -  -B-'\_A(ed  -  Kjfi  -  ed)  -  K„(B  -  ed))  -f]  (3  - 19) 

The  control  angle  (<5)  is  a  function  of  the  pitch  angle  (6)  ,  desired  pitch  angle 
( dd ),  angular  velocity  ( 6 ),  desired  angular  velocity  (6d),  desired  angular  acceler¬ 
ation  (6d),  the  position  feedback  gain  Kp  ,  the  velocity  feedback  gain  Kv,  and  the 
matrices  i.e.,  A,  B,  F.  Since  the  thrust  force  T  is  included  in  F  and  B>  the  mag¬ 
nitude  of  the  thrust  force  will  thus  affect  the  control  angle  S  .  Kv  and  Kp  are  ad¬ 
justable  to  obtain  the  desired  response  of  the  pitch  angle  6. 
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IV.  COMPUTER  SIMULATION 


A.  SIMULATION  OBJECTIVES 

In  literature,  the  application  of  forces  nas  been  limited  to  gravity.,  aerodyna¬ 
mic  forces  and  thrust  [Ref.  2].  No  damping  has  been  implemented.  The  thrust 
and  aerodynamic  forces  can  be  found  fairly  accurately  with  standard  test  and 
design  procedures  such  as  static  firings  to  obtain  thrust  versus  time  profiles  for 
the  engine  and  wind  tunnel  measurements  to  determine  the  aerodynamic  forces. 
Gravitational  forces  can  be  calculated  from  the  knowledge  of  the  missile's  posi¬ 
tion  relative  to  the  earth.  The  mass  including  fuel  can  be  estimated  from  know¬ 
ledge  of  the  missile's  weight  before  and  after  burnout  (from  the  measurements), 
and  by  using  a  mathematical  relationship  (often  linear)  for  the  decrease  in  missile 
mass  over  the  engine  burntime. 

In  this  work,  the  missile's  weight  is  assumed  constant  and  the  aerodynamic 
forces  are  zero  The  deformations  resulting  from  the  structural  flexibility  have 
been  assumed  small  and  small  control  angle  assumptions  are  used  in  the  rigid- 
body  controller  design. 

The  primary  purpose  of  this  study  is  to  complete  the  simulation  of  the  flexible 
missile  in  2-D  motion,  where  the  bending  effect  is  considered  to  be  the  only  flex¬ 
ibility  source  and  a  rigid-body  controller  as  developed  in  the  last  section  is  in¬ 
cluded. 

B.  SOLUTION  TECHNIQUE 

Many  numerical  integration  techniques  can  oe  applied  in  the  solution  of  the 
equations  of  the  motion.  The  main  consideration  in  the  selection  of  the  inte¬ 
gration  technique  is  the  size  of  the  time  step  necessary  to  integrate  the  equations 
of  motions  with  numerical  accuracy  and  stability. 

The  type  of  the  equations  of  motion  (2-20,  2-21)  in  this  research  permits  the 
application  of  the  Sequential  Integration  Method  [Ref.  6].  The  linear  equations 
of  small  motions  are  integrated  implicitly  and  the  large  motion  solutions  are  ob¬ 
tained  using  explicit  integration  method.  The  implicit  methods  are  effective  for 
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linear  systems  with  high  frequencies  and  the  explicit  methods  are  effective  for 
solving  nonlinear  systems  with  low  frequencies.  The  implicit  method  is  especially 
effective  for  linear  systems  having  a  wide  range  of  frequencies  of  which  only  the 
lower  frequencies  are  excited.  The  size  of  the  time  step  need  only  be  chosen  to 
make  the  solution  of  the  excited  modes  sufficiently  accurate.  Explicit  methods 
are  effective  for  large  scale  systems  with  low  frequencies.  Because  of  the  low  fre¬ 
quencies,  the  size  of  the  time  step  that  we  choose  iur  stability  need  not  be  small. 
In  addition,  the  explicit  method  does  not  need  iterative  procedures,  which  are 
time  consuming  for  non-linear  systems. 

C.  THE  COMPUTER  SIMULATION  CODE 

A  high  level  computer  language,  i.e.,  MATLAB,  was  chosen  to  simulate  the 
missile  system.  The  MATLAB  was  designed  for  matrix  operations.  The  simu¬ 
lation  code  was  developed  with  modular  MATLAB  routines. 

The  simulation  code  can  be  divided  into  three  levels  :  LEVEL  1  (an  overview) 
separates  the  code  into  primary  portions  of  INITIALIZATION,  PLANT  DE¬ 
SCRIPTION,  INTEGRATION,  SYSTEM  CONTROL  and  OUTPUT.  LEVEL 
2  facilitates  several  subroutines  for  LEVEL  1.  A  listing  of  MATLAB  routines 
required  for  manipulation  in  LEVEL  2  is  placed  in  LEVEL  3. 

D.  SIMULATION  RESULTS 

The  computer  simulation  was  performed  with  variable  parameters  which  de¬ 
termine  missile  speed  and  controller  bandwidth.  These  parameters  include  force 
T,  A„  and  Kf. 

The  simulation  outputs  will  be  presented  in  large  motions,  small  motions  and 
control  angle.  The  large  motions  are  X0,  J'0  and  and  the  small  motions  are  v(0) 
and  <3>(0).  The  initial  condition  for  all  runs  was  the  pitch  angle  of  45°.  The  sim¬ 
ulation  work  was  divided  into  two  areas.  First  a  simulation  was  performed  for 
the  rigid  missile  using  the  rigid-body  controiler.  These  results  were  used  as  a 
baseline  for  comparison  with  the  results  of  the  flexible  missile  using  the  rigid -body 
controller  system.  Second  a  simulation  was  performed  for  the  fiexible  missile  and 
rigid-body  controller  system.  Only  the  trajectory  control  will  be  discussed  in  the 
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analyze  of  the  dynamics  of  the  flexible  missile.  The  desired  trajectory  was  speci¬ 
fied  as 


45’ 


6  “ 


45°  -  112.5/ 


when  0.1  £  t  <  0.1 
when  0.1  £  t<  0.4 
when  t  £  0.4 


The  control  angle  is  assumed  to  be  limited  to  +/-10  degrees  (small  angle). 
This  restriction  puts  a  saturation  line  to  the  control  input. 

Fig.5  represents  the  graphical  results  of  the  desired  and  actual  trajectory  tor 
the  rigid  missile  an<l  rigid-body  controller.:  The  force  (T)  and  controller  band¬ 
width  (cuj  were  30000  N  and  3  ra d/'s  respectively.  The  dark  black  line  and 
dashed  line  represent  the  desired  and  actual  trajectory,  respectively.  The  control 
angle  (<5)  is  shown  in  Fig  .6  and  Figures  7-8  present  the  base  positions  (A^F0)  in 
large  motion. 

After  presenting  results  for  the  rigid  missile  with  rigid-body  controller,  the 
dynamic  behavior  of  the  flexible  missile  with  rigid- body  controller  will  be  studied 
next.  The  force  (T)  and  controller  bandwidth  (w*)  again  were  30000  N  and  3 
rad /$,  respectively.  Fig-9  shows  the  desired  and  actual  trajectory.  The  control 
angle  Is  presented  In  Fig.10.  After  0.1  sec,  the  missile  Initially  needs  a  control 
angle  which  is  less  than  10  degrees.  The  effects  of  small  motion  can  be  seen  on 
control  angle  clearly.  The  base  deflection  and  slope  are  shown  in  Figures  11-12. 
There  is  no  small  motion  effects  on  flexible  missile  between  0-0.1  sec  because  of 
zero  control  angle  and  ho  force  components.  After  0.1  sec,  the  small  motions  are 
excited  and  have  an  amplitude  of  10-*.  Figures  13-14  present  the  base  positions 
of  the  flexible  missile  In  large  motion.  The  large  motion  behaves  like  the  rigid- 
body  motion,  which  implies  that  the  coupling  effect  between  the  large  motion  and 
small  motion  is  small.  The  dynamic  behavior  of  the  flexible  missile  is  dominant 
by  the  large  motion  in  this  case  since'  the  bandwidth  of  the  controller  is  low  and 
small  motion  is  not  significantly  excited. 


RIGID  MISSILE  CONTROL  SYSTEM 


Figure  5.  The  desired  and  actual  trajectory  for  the  rigid  missile  with  rigid-body 
controller  (T  •  30000  N,  <on  -  3  rad/s) 


TIME(S) 


The  control  angle  for  the  rigid  missile  with  rigid-body  controller  (T 

30000  N,  ct)„  —  3  »ad/s) 
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Figure  7.  The  large  motion  position  X  for  the  rigid  missile  with  rigid-body  controller 
(T  -  30000  N,  to„  ™  3  rad/s) 
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Figure  8.  The  targe  motion  position  Y  for  the  rigid  niissite  with  rigid-body  controller 
(T  -  30000  N,  -  3  rad/s) 
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Figure  9.  The  desired  and  actual  trajectory  Tor  the  flexible  missile  with  rigid-body 
controller  (T  *=  30000  N,  w,  ™  3  rad/s) 
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f  igure  10.  The  control  angle  for  the  flexible  missile  with  rigid-body  controller  (T 

-  30000  N,  to,  -  3  rad/s) 
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Figure  1!.  The  small  motion  position  v(0)  for  the  flexible  missile  with  rigid-body 
controller  (T  “  30000  N,  a>„  ■=  3  rad/s) 
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Figure  12.  The  small  motion  position  phi(O)  Tor  the  flexible  missile  nith  rigid-body 
controller  (T  —  30000  N,  a>„  =  3  rad/s) 
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Figure  14.  The  large  motion  position  Y  tor  the  flexible  missile  nith  rigid-body  con¬ 
troller  (T  «  30000  N,  (o„  -  3  rad/s) 
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Two  kinds  of  tests  were  performed  to  further  the  study  of  the  dynamic  be¬ 
havior  of  the  flexible  missile. 

The  first  test  was  to  study  dynamics  due  to  the  increase  of  the  control  system 
bandwidth  (con).  The  desired  and  actual  trajectory  of  the  flexible  missile  with 
rigid-body  controller  without  saturation  on  the  control  angle  is  presented  in 
Fig. 1 5.  The  force  (T)  is  still  30000  N  while  the  controller  bandwidth  (co„)  is  in¬ 
creased  to  300  rad/s.  As  the  controller  bandwidth  is  increased,  the  gap  between 
controller  bandwidth  and  system  natural  frequency  is  smaller.  The  system  will 
be  unstable  at  controller  bandwidths  close  to  the  system  fundamental  frequency 
which  can  be  as  high  as  270  rad/s  in  this  case.  Note  that  the  fundamental  fre¬ 
quency  of  the  simple  beam  with  free-free  boundary  conditions  wras  calculated  as 
coM  —209.2053  rad/s,  and  the  missile's  fundamental  frequency  will  be  increased 
due  to  the  coupling  between  large  and  small  motion.  The  control  dynamics  in¬ 
terfere  with  the  structural  dynamics.  The  higher  control  frequency  with  the 
bandwidth  of  300  rad/s  significantly  excites  the  small  motion  of  the  structure. 
The  graphs  of  control  angle  and  the  small  motions  and  large  motions  of  the  mis¬ 
sile  base  (Figures  1 6- 1 7-  i  8- 1 9-20)  show  the  unstable  state. 

The  bandwidth  of  the  controller  must  be  set  much  lower  than  the  funda¬ 
mental  frequency  of  the  missile  in  order  to  use  the  rigid-body  control.  This  implies 
that  the  pitch-angle  response  speed,  which  is  determined  by  the  control  band¬ 
width  is  limited.  To  keep  the  response  speed,  the  missile  structure  must  be  de¬ 
signed  sufficiently  rigid  to  possess  a  high  fundamental  frequency.  On  the  other 
hand,  the  missile  payload  will  affect  the  natural  frequency  of  the  missile  structure 
with  the  heavier  the  payload,  the  lower  the  fundamental  frequency  of  the  missile. 
Therefore,  the  payload  must  oe  limited  to  achieve  high-speed  control  response. 
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Figure  15.  The  desired  and  actual  trajectory  for  the  flexible  missile  uith  rigid-body 


controller  without  saturation  on  the  control  angle  (T  -  30000  . 1,  gj„ 


*  300  rad/s) 
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I  igure  16.  Tlw  control  angle  for  the  flexible  missile  with  rigid-body  controller 
without  saturation  on  the  control  angle  (T  -  30000  N,  c on  -  300 
rad/s) 
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I  igure  18.  The  small  motion  position  phi(O)  for  the  flexible  missile  rigid-body  con¬ 
troller  without  saturation  on  the  control  angle  (T  *  30000  N,  a>„ 
300  rad/s) 
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I  isure  19.  The  large  motion  position  X  for  the  flexible  missile  rigid-body  controller 
without  saturation  on  the  control  angle  (T  •=  30000  N,  co„  »  300 
rad/s) 
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Figure  20.  The  large  motion  position  Y  for  the  flexible  missile  rigid-body  controller 
without  saturation  on  the  control  angle  (T  -  30000  N,  co  «  300 
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Fig-2 1  presents  the  desired  and  actual  trajectory  of  the  flexible  missile  using 
rigid-body  controller  with  a  saturation  (  +  /-10  Degrees)  on  the  control  angle.  The 
force  (T)  and  controller  bandwidth  (co„)  are  30000  N  and  200  rad/s,  respectively. 
A  better  trajectory  tracking  is  obtained,  however  the  actual  trajectory  is  oscillat¬ 
ing  about  the  desired  trajectory  which  is  mainly  due  to  the  switching  of  the  con¬ 
trol  angle  between  the  saturation  lines.  Fig. 22  shows  the  control  angle,  and 
Figures  23-24  show  the  small  motion  displacement  and  slope  of  the  base  of  the 
flexible  missile.  It  was  observed  that  the  switching  of  the  control  angle  between 
+  /-10  degrees  causes  increased  excitation  of  small  motions.  Figures  25-26  pre¬ 
sents  the  position  of  missile  base  in  large  motion  where  the  large  motion  is  no 
longer  dominated  by  the  rigid-body  motion. 

The  second  test  was  to  study  the  dynamics  of  the  flexible  missile  due  to  an 
increase  of  the  missile  speed.  To  increase  the  speed  of  the  missile,  the  thrust  force 
was  increased.  Fig. 27  presents  the  actual  and  desired  trajectory  of  flexible  mis¬ 
sile.  The  force  was  60000  N  and  controller  bandwidth  was  3  rad/s.  As  seen  from 
Fig. 27,  the  pitch-angle  response  is  determined  directly  from  the  bandwidth  of  the 
controller.  The  increased  missile  speed  does  not  change  the  pattern  of  the  pitch- 
angle  response.  The  control  angle  of  the  flexible  missile  is  shown  in  Fig. 28.  The 
control  angle  is  affected  by  the  increased  speed.  When  the  speed  is  increased,  the 
control  angle  gets  smaller  because  a  smaller  control  angle  generates  a  sufficient 
correction  to  the  pitch  angle.  The  increased  missile  speed  has  little  effect  on  small 
motions  (Figures  29-30).  Figures  31-32  show  the  large  motions. 
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MISSILE  CONTROL  SYSTEM 


(SHDa)NOUISOd  HVinONV 


Figure  21.  The  desired  and  actual  trajectory  for  the  flexible  missile  using  rigid-body 
controller  with  saturation  on  the  control  angle  (T  *■  30000  N,  co„  — 
200  rad/s) 
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TIME(S) 


MISSILE  CONTROL  SYSTEM 


figure  22.  The  control  angle  for  the  flexible  missile  using  rigid- body  controller  with 
saturation  on  the  control  angle  (T  *“  30000  N,  (x>„  "  200  rad/s) 


38 


TIME(S) 


(W)NOIJLISOd  INHWHOVldSICI  NOHOW  HVWS 

gure  23.  The  small  motion  position  v(0)  for  the  flexible  missile  using  rigid-body 
controller  with  saturation  on  the  control  angle  (T  -  30000  N,  u)n  = 
200  rad/s) 
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MISSILE  CONTROL  SYSTEM 


1  igure  2u.  The  (urge  motion  position  X  for  the  flexible  m  issile  using  rigid-body 
controdet  with  saturation  on  the  control  angle  (T  =*  30000  N,  tx>„  — 
200  rad/s) 
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MISSILE  CONTROL  SYSTEM 


(W)NOllISOd  A  TVBOIO 

Figure  26.  Tlie  large  motion  position  Y  for  tiie  flexible  missile  using  rigid-body 
controller  with  saturation  on  the  control  angle  (T  =  30000  N,  o)n 
-  200  rad/s) 
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Figure  27.  The  desired  and  actual  trajectory  for  the  flexible  missile  using  rigid-body 
controller  with  Increased  speed  (T  —  60000  N,  con  “  3  rad/s) 
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MISSILE  CONTROL  SYSTEM 


Figure  28.  The  control  angle  for  the  flexible  missile  using  rigid-body  controller  with 
increased  speed  (T  =*  60000  N,  a>„  =*  3  rad/s) 


TIME(S) 


(W)NOIJJSOd  ±N3W33Y33Sia  NOIJLOIM  T1VVNS 

The  small  motion  position  v(0)  for  the  flex'. ble  missile  using  rigid-body 
controller  n  ills  increased  speed,  (T  «  uG'jOO  N,  to,  —  3  rad/s) 


xi(H 


(w)nouiso<j  adds  nouow  tiviais 


m  ""  !“sl11”'  nm  'or  the  flexible  uslllg  rigW  ~ 

conltoher  nitli  Increased  speed  (T  »  6QOOO  N,  u)H  -  3  rad/s) 
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MISSILE  CONTROL  SYSTEM 


Figure  3i.  The  large  motion  position  X  for  the  flexible  missile  using  rigid-body 


controller  with  increased  speed  (T  =  60000  N,  a)„  —  3  rad/s) 
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TLME(S) 


Figure  31  The  large  motion  position  Y  fo^  the  flexible  miss^e  using  rigid-body 
ccntroller  with  increased  speed  (7  FGCO0  N,  co„  “  3  rad/s) 
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V.  SUMMARY 


A.  CONCLUSIONS 

The  development  and  computer  simulation  is  presented  for  a  bending  flexible 
missile  with  a  rigid-body  controller  system  moving  in  a  2-D  coordinate  frame.  In 
this  /research,  the  dynamic  model  , has  been  developed  through  the  Equivalent 
Rigid  Link  System^  utilizing  Lagrangian  dynamics  to  obtain  a  type  of  system 
equations  of  motion  suited  for,  Sequential  Ktiegration  Method  that  integrates 
large  motion  explicitly  and  small  motion  implicitly.  The  spatial  finite  element 
discretization  of  missile  structure  and  the  application  of  truncated  natural  modal 
responses  provide  an  approximate  solution. 

The  analysis  and  simulation  were  performed  to  understand  the  dynamic  be- 

«» 

havior  of  a  flexible  missile  using  a  rigid-body  controller.  It  was  found  that  the 
controller  bandwidth  must  be  much  lower  than  the  fundamental  frequency  of  the 
missile  in  order  to  use  the  rigid-body  controller.  The  payload  will  affect  the  na¬ 
tural  frequency  of  the  missile  structure  i.e,  when  the  payload  is  incieased,  the 
system  fundamental  frequency  will  be  decreased.  The  payload  must  then  be  lim¬ 
ited  to  achieve  high-speed  response.  In  order  tc  increase  the  payload  and  main¬ 
tain  high-speed  control  response,  a  flexible-body  controller  is  needed.  /’ 

0.  RECOMMENDATIONS 

Areas  which  remain  to  be  investigated  include  : 

1. Add  the  payload  and  aerodynamic  effects  to  the  model. 

2.  Design  and  study  the  dynamic  behavior  of  a  flexible  missile  with  a 
flexible-body  controlier  in  2-D  motion. 

3.  build  a  flexible  m'ssile  in  a  laboratory  scaie  and  obtain  experimental  ata. 

4.  Design  and  rimulate  a  control  system  for  flexible  missiles  in  3-D  motion. 
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APPENDIX  A.  DERIVATION  OF  THE  MODE  SHAPE  RESPONSE 
MATRIX  FOR  THE  FLEXIBLE  MISSILE 


In  this  study,  the  natural  mode  shape  functions  of  a  beam  are  used  to  repre¬ 
sent  the  flexural  motion  of  the  flexible  missile.  Only  the  first  two  mode  shapes  are 
used.  The  flexible  missile  is  modeled  as  a  continuous  Euler-Bernoulli  free-free 
beam,  neglecting  shear  deformation  and  rotary  inertia  effects. 

The  bar  (Fig.33)  has  system  parameters: 


Figure  33.  The  Bar  in  Flexure 


v(x,t)  “Transverse  displacement  at  any  point  x  and  time  t 

f(x,t) -Transverse  force  per  unit  length 

m(x)“The  mass  per  unit  length 

El(x)  “Flexural  rigidity 

E  '  Young's  modulus  of  elasticity 

I(x)  “The  cross-sectional  area  moment  of  inertia 

Q(x,t) “Shearing  rorce 

M(x,t)  “  Bending  moment 


Fig.34  shows  the  free  body  diagram  corresponding  to  a  bar  element  of  length 


Figure  34.  Free  Body  Diagram  Corresponding  to  a  Bar  Element 


The  force  balance  of  the  free  body  is 


CGC*,0  +  —  <it]  -  Q(x,t )  +J{xft)dx  -  m{x)dx  - - 

dx 


(A-  1) 


The  moment  equation  of  motion  about  the  axis  normal  to  x  and  y,  ignoring 
the  inertia  torque  associated  with  rotation, 


5! 


lM(x,t)  + 


dM(x,t)_  _  M(x,t)  +  [Q{x,i)  +  — dx]dx 


dx 


8x 


+  J{xd)dx- y-  =  0 


(A -2) 


Canceling  appropriate  terms,  ignoring  terms  involving  second  order  terms  in 
dx  and  combining  Eqs.  (A-l)  and  (A-2), 


d2M(x,t) 


dx 


+Axd)  -  M(xJ) 


cy{*d) 


j  ■  ^r»v  v*-/  2 

2  dr 


(A-2) 


Eq.  (A-3)  relates  the  bending  moment  M(x,t),  transverse  force  f(x,t)  and 
bending  displacement  y(x,t). 

The  relation  between  the  bending  moment  and  the  bending  deformation  is 


M(x,t)  =  £I(x,t ) 


d2y(x,t) 

dx2 


(A-  4) 


Inserting  Eq.  (A-3)  into  Eq.  (A-4),  we  obtain  the  differential  equation  for  the 
flexural  vibration  of  a  bar, 


C?*2 


(£/(*) 


dx 


)  +/M  =  A/(jc,r) 


d2y{x,t) 


2  >  ■  j  _  ? 

dt 


(A  -  5) 


The  bending  moment  and  shearing  force  of  the  free  ends  (x~0,  x  =  L)  are 
zero. 


£/(*) 


d2y{xj) 


dx 


2  'x=0 


Uo  -  0 


„„  ^ d  y(xA)  , 

£/(x)  — --  \X=L  =  0 
dx 


(A-  6) 


J^lei{x)1^!L^  =  0 


d  r  r'  n  \  d2y(xd) 

—  [£/(*)  -p- ]  Uz.  -  0 


{A-l) 


Eqs.  (A-6)  and  (A-7)  are  called  natural  boundary  conditions. 
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Considering  the  free  vibration  characterized  by  f(x,t)=0,  using  separation  of 
variables  method  and  simplifying  Eq.  (A-5), 


dx2 


dxi 


]  =  w*m(x)  Y{x) 


(A-  8) 


Simplifying  the  Eq,  (A-8),  for  EI(x)  -  constant 

d*Y(x)  .  ,  x 
- Y~~fiiY(x)  =  Q 


dx' 


where  /?4  = 


w2w(x) 

Efixf 


(A-  9) 


The  boundary  conditions  require  that, 
At  x  =  0  (base) 

Mi  _0 

Uo”° 


<rr{x) 


dx 


3  lx=0 


=n=0 


{A  ~  10) 


At  x  =  L  (tip) 


d2  Y(x) 


dx 


2  lx =L 


=  0 


d‘Y(x) 
dx 3 


=  0 


(-1-11) 


The  general  solution  of  Eq.  tA-9), 

F(.t)  -  C\  sin  {fix)  4-  C2  cos(Px)  +  C3  sinh  (Px)  +  C'4  cosh(/Jx)  (A  -  12) 

Taking  the  derivatives  of  Eq.  (A- 12),  and  substuting  boundary  conditions, 
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K(x)  =  C](  sin  (lx  4-  sinh  (lx )  4-  C2(  cos  Px  4-  cosh  px) 


-  13) 


Solving  Eq.  (A- 13)  yields, 

cos  pL  cosh  pL  =  1  (A  —  14) 

The  first  five  consecutive  roots  of  this  equation  are; 

fit  =  0.0 

p2  =  4.712388 

/?,  =  7,853981 

/?4=  10.995574 

/?,=  14.137166 

Eq.  (A- 13)  is  now  written  in  the  following  form, 

Yr(x)  =  Cr(  cos  pyX  4-  cosh  PyX)  4-  ( sin  PyX  4-  sinh  pyX)  r—i,2  (4-15) 


where: 

^  sin  PrL  -  sinh  (lrL 
r~  -  COS  PyL  Jr  cosh  pyL 


(A  -  16) 


The  transverse  displacement  v(x)  and  slope  <l>(x)  can  be  represented  in  the 
following  forms  respectively, 

v(x)«  iarYr(x)  (A  -11) 

r— 1 

cP(x)  -  — -  (A  —  18) 

ox 


v(jr)  =  ax(Ci(  cos  PiX  4-  cosh  P .x)  4-  ( sm  fixx  4-  sinh  pxx)) 

+  a2 (C2(  cos  p2x  +  ccsh  p2x)  +  (  sin  P2X  +  sinh  p2x))  ( A  —  1 9) 

<D(x)  -  a^CiPi x  +  sinh  Pxx)  4  /?,(  cos  fixx  4-  cosh  Pxx)) 

4-  a2(C2p2(  -  sin  p2x  +  sinh  P2x)  -r  ( cos  p2x  4-  cosh  P2x))  {A  —  20) 
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Substituting  the  boundary  conditions  into  shape  function  gives, 


v(0)  —  2fljCj  3*  '2.cl-^C2 

(T  —21) 

<I>(0)  =  2«sr1^|  4-  2i72/?2 

(/I  -  22) 

The  modal  amplitudes  at  and  a2; 

*>-<2, :c,  +  c:5)vl0)+c,£  <t(0) 

-  23) 

=  ~  l,(0)  +  T <!>(0) 

(/l  -  24) 

v  =  tfv(O)  4  6<t>(0) 

(4  -  25) 

a  =  F1r,p)4f3K2W 

-  26) 

*-*i>i(*)  + 

(A  -  27) 

F  1  l  C 2^1 

f'~2c;+  c>e 

(A  -  28) 

Cic 

(.4  -  29) 

F-  J' 

(A  -  30) 

f^~~e 

(4  -  31) 

C2P1 

E  —  2P2  —  2  -- “~- 

(A  -  32) 

Substituting  Eqs.  (A-15),  (A-16),  (A-26),  (A-27),  (A-28),  (A-29),  (A-30), 
(A-3I),  (A-32)  into  Eq.  (A-25)  yields  an  expression  for  the  transverse  displace- 


ment  of  a  flexible  missile  as  a  function  of  the  missile  base  nodal  displacements, 
v(0)  and  <I>(0), 

v(x,r)  —  [/^(C^  cos  Pyx  4-  cosh  pyx)  4-  ( sin  p{x  4-  sinh  p{x)) 

4  F3(C2(  cos  p2x  4  cosh  p2x)  4  (  sin  p2x  4  sinh  P2x))]v(0) 

4  [F2(Ci(  cos  Pyx  -f  cosh  pyx)  4-  ( sin  pxx  4  sinh  Pyx)) 


4-  F4(C2(  cos  P2x  -i-  cosh  P2x )  -1-  ( sin  p2x  4-  sinh  -  33) 


This  expression  is  differentiated  twice  to  obtain  v"(„v),  which  is  necessary  for 
the  calculation  of  the  potential  energy  due  to  deformation  and  theoretical  strain. 

v"(x,r)  =  (fy(Cyp\(  -  cos  Pxx  4*  cosh  /?,)  4-  p\{  —  sin  Pyx  4-  sinh  /?,x)) 

4-  F3{C2pl(  -  cos  P2x  4  cosh  P2)  4-  p\{  —  sin  p2x  4-  sinh  p2x)))v( 0) 

4-  (F2(Clp‘]{  -  cos  pyx  4-  cosh  Py)  4-  p\{  -  sin  pxx  4-  sinh  pyx)) 

4-  F4{C2Pl  —  cos  p2x  4-  cosh  P2)  4-  P2(  -  sin  P2x  4-  sinh  p2x)))<t>{ 0)  (/I  —  35) 


Now  substitution  of  v(x)  into  the  3x1  deformation  vector,  d  ,  yields  the  3x2 
shape  i'un;tion  matrix,  </>,  and  the  2x1  nodal  displacement  vector,  V , 


d  =  <pv  • 


0 

' o  o' 

r*  ,  i 

0 

0  0 

|  vp) 

— 

[<0(0) 

v(x)_ 

g  k 

L.  -J 

( A  -  36) 


The  shape  function  matrix  is  now  in  a  form  convenient  for  computer  coding. 
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APPENDIX  B.  DERIVATION  OF  THE  EQUATIONS  OF  MOTION  FOR 

THE  FLEXIBLE  MISSILE 

The  ERLS  is  used  to  separate  the  flexible  missile's  motion  into  a  small  motion 
and  a  large  motion.  The  large  motion  generalized  coordinates  are  defined  as  X0, 
K0,  theta  (0)and  the  small  motion  generalized  coordinates  as  U.  The  Lagrangian 
Dynamics  are  very  helpful  in  the  derivation  of  the  equations  of  motion  of  com¬ 
plex  systems.  U  represents  a  column  vector  which  has  two  elements.  First  one  is 
v(0)  which  is  the  nodal  displacement  of  the  flexible  missile's  base  and  d>(0)  which 
is  the  slope  cf  the  flexible  missile's  base.  We  can  represent  the  vector  of  the  gen¬ 
eralized  coordinates  by  Eq.  (B-l), 

?-[  *0  n  e  ]r  £  =  0(0)  0(0)  ]r  (5-1) 

These  generalized  coordinates  will  be  used  in  the  application  of  the  Lagrange 
equations  to  the  equations  of  motion  for  both  small  and  large  motion  coordinate 
stems.  The  Lagrangian  equations  are  defined  in  Eqs.  (B-2)  and  (B-3). 


d  r  dKE  dKE  ,  dPE  r- 


(/ =  1,2,3) 


(B  -  2) 


d  r  dKE 
dt  l 

dU 


dKE  , 

- Z1 - r 

dU 


3  PE 
dU 


(B  -  3) 


where  is  a  component  of  £ 

Using  ERLS  we  can  define  the  global  position  of  the  base  position  in  terms 
of  the  local  position  vector  (r),  a  deformation  vector  ( d ).  and  a  coordinate  trans¬ 
formation  matrix  (W)  with  Eq.  (B-4), 

R  =  B'(r  +  d)  ( B  -  4) 
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Continuing  in  the  development  of  the  equations  of  motions,  we  determine  the 
derivative  with  respect  to  time  for  global  positions  in  order  to  obtain  the  kinetic 
energy  expressions. 


R~W(?  +  d)  +  Wd  (B  -  5) 

The  kinetic  and  potential  energy  expresions  follow.  In  this  development  we 
have  separated  the  large  and  small  motion  energy  contributions  and  will  present 
them  seperately. 


KE  -  ~ \RrRdm 

PE  =  ~ \UTVTCVUdx  -  \Rrgdm 


{B-  6) 
(5-7) 


Redefining  d  in  terms  of  the  shape  matrix  <p  (derivation  of  d  is  presented  in 
Appendix  A  ),  we  can  rewrite  Eq.  (B-6)  as 


KE  =  j-  J([Pk(r  4  d)  +  WdflWU  4  d  +  l,  d])dm  ( B  -  8) 

KE  =  -^\{7TWTWf  +  2rTWTW(pU  +  2  TV/T\V<pU 

+  UT(j)TlVTlV(t>U  +  2 UT<pTWTW4,U  +  UT<t>TWT\V4>UT)drn  ( B  -  9) 

Continuing  in  the  development  of  the  equations  of  motion,  we  first  express 
the  derivative  of  the  kinetic  energy  with  respect  to  the  time  rate  of  change  of  the 
large  motion  coordinate  £,  Eq.  (T?-13),  and  then  determine  the  time  rate  of  change 
of  this  expression,  Eq.  (B-14;.  Since 
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(5  —  10) 


■*o“ 

V 

^0 

11 

^3 

IV  2 - 


n\  = 


dx 

div 

dY 

div 

de 


dw 


Wtl  =  ~  =  Partial  derivative  w.r.t.  X 


8VV 

3*2 

d\V 


=  Partial  derivative  w.r.t.  Y 


=  Partial  derivative  w.r.t.  0 


IK,.  =  Partial  derivative  of  W  w.r.t.  X,  Y,  6  i  =  1,2,3 


dKE  +  +  Sly 


i  ,T 


K  2  ' '  K 


Ki 


w<i>u 


+l7T\VT——  4>U  4-  l7T  ( lK</>£/  4- 


(?S| 


l'V^)U 


)T  lT  li/T  dlK 


:rr,7  dll 


,7" 


4-  U‘  <p'  IV '  4>U  +  21  J 1  4> 1  -^r  ~  W<t>U  )  dm 

dii  dti 


Putting  Eq.  (B-12)  into  Eq.  (B-ll)  and  simplifying 


d\V  d\V 


d'Zt  dt> 


d\VT  _  dlV7  _  ,,/r 
- 


'w 


dX£ 


=  j (  r  rW',J  1P7  4-  \v]w<pU  f  rTWT\V^U 


+  rT wj W4>ij  +  UT<pT w] W<pu  4-  VT4>T iv J IVtfiU  )  dm 


(B  —  11) 


(5-12) 


(5-13) 


-jL  ( ML  )  =  j(  ytwliv?  +  rrivj iW  4-  rriv;3iv<t>v 


+  rT Wj lV<pU  4-  P'WJ W(f>U  4-  rT WT^<f>U 

+  7TWTW^V  4-  fTWTW^V  +  7TIV~3IV(PU 
4-  rTW J V/<f>U  4-  rT w] W6U  4-  UT(f>  Twj IV$U 

4-  WtvfyvfV  4-  ur(f>TivJiv<f>u  4-  uT4>TivJiv<t>u 
4-  Ur4>TlvJlV<t>U  4-  uWlVjlVtU  4-  UT<t>rtvflV<l>U 

4-  4/ V  IFJ (*>£/ )  am  {3  ~  14) 

Finally  rve  express  the  derivative  of  the  kinetic  energy  with  respect  to  the 
laige  motion  position  as 

-----  =  f(  7TWT,\Vr  4-  7TW1iW$d  -r  TtiVtW}4>V 
4-  rr  1^3 «'$>{/  4-r7W,7V^i/  4-  UT4>T\vllV<t>U 

4-  UT<i>TWT3W<t>U  4-  LJ1 4>r\VTWJ4>U  4-  £/r<prtt'?  )  dm  ( B  ~  1 5) 
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(. B  -  16) 


A,  (ML  |  „  ML  =  j(rVJif7  +  7TwJw<pu  +  2rr^J^c? 


+  rrf^T^Z/  +  rrlvjlV4>U  +  UT<t>TwJW(f>U 


H  2LtT4>TwJw$U  4-  UT<f>TwJlV<f>U )  dm 


Completing  our  development,  the  expressions  for  the  potential  energy  are 
presented  organized  similarly  to  the  previous  material, 


=-|(Tf  <f>U)rwJgdm  =  -  jP'iVjgdm  -  j  UT<pTWTigdm  (B-  IT) 


The  small  motion  equations  are  handled  in  a  similar  way,  and  are  easy  to 
derive  than  large  motion. 


^  j(4>TlVTl V?  f  4>tWtW4>U  4-  <t>TWTW<}>U)dn 7  (3-  IS) 

dUT 


ML  =  j  (<f)TlVTlv7  +  4>t\ VTW<t>U  +  4>TWTW<l>U)dm  (B  —  19) 

QL‘ 


~dt  ~ ) =  WTWr  +  <}>TWTW T  +  <}>TWTW<j>U  4-  <\>TwTw<t>u 
dVT 


- 'r2<bTWTW<t>U  +  <}>TWTW4>U  +  4>tWtW4>U  )  dm 


A 

dt 


ML) 

dUT 


-  =  \(<t>TWTWr  4-  <f)TWTiV<t>U 

duT 


(£  -  20) 
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{B  -  21) 


H  2 <prV/TW<t>d  4  4>TWfW<i>U  )  dm 
Potential  energy, 

~~  =  frrcrt;<i>:-  -  rfm  (i?  -  22) 

o'/7 

We  can  simplify  these  expressions  further  by  (he  substution  of  the  second 
time  derivative  of  the  coordinate  transformation  matrix  W.  The  contributions  to 
the  coriolis  and  centripetal  forces  represented  by  \VT  ,  termed  the  residual  accel¬ 
erations.  The  contribution  to  the  general  forces  is  represented  by  Wti . 

W  =  Wr  4  I  W£t  (B  -  23' 

/=i  ’ 

For  large  motion, 


f5M<  L  wfi,  +  w,)7  +  ?rwJ{  i  wj;,  +  wr)t,a  +  irT ivjv^u 
J= 1  J" 1 


3  _  _  _u  3 

+  rr(  I  1F>.  -4-  ~yVr)Wj<f)U  +  rrH'jWpU  +  Ur4>rivj(  £  +  1F,)<£(/ 

>=i  1  '  ’  ;=i 


+  ;.UT<t>TlvjlV4>U  +  Uj<hTw}  W<t>U  -  r7' wjg  ~gTVV^U)dm  ■=  (J9  -  24) 

For  small  motion, 

\{$TWT{Z  WjZj  4  lFr)7  +  ^rlFr(|  +  JF,)<££/  +  2r<pr WTW<pU 

j— l  y=i 


62 


(B  -  25) 


+  4>TwTv/fu  +■  rrcru  -  <bTwT%  -  ?u 

Collecting  tenns  <md  arranging  the  coefficients  the  two  Lagrange  equations 
can  bt  written  as  the  equations  of  motion  for  the  large  and  small  generalized  co¬ 
ordinates.  It  is  these  equations  which  must  be  solved  by  computer  simulation 
code. 

For  the  large  motion, 


3 

I 

;= i 


F'Wj  W/dm  + 
frrlvJWj<t>Udm  + 
\UT$Tl"JWj4>Udm  + 
jUT</)T IV J IV Jr  dm 


\rTWT:  W(j)dm  + 

\UT<t>'rwJw<t>'im 


U  = 


FfJ  -  \7rwrt  Wfdm  -  \7T wj W,4>Udm  ~  2 \7T w] WQUdtn 
jUT(j)TwJ Wr7dm  -  jUT4>Twj W4Vdm  ~2j  UT<j>TwJlV<pUdm 

4-  \rT IV Jg dm  -V  jgTW j(j)Udm 


(B  -  26) 


For  the  small  motion, 

3 

X  [\4>TlVTW/dm  +  \<t>TWTWJ4)Udm\ij  +  [\<i>TWTW4>dnvU  4- 

y=i 

[f2  $TWrW4>dm\U  +  [j(pTWr<t>dm  4-  jl^CTWJt/  = 

|  <j)  TWTg  dm-j<j>TW7Wsdm  +  Fu  (5~27) 
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By  cefming 


Mqq(iJ)  =  j 7TW*W~rdrr  4-  \7Tlvjwj(pUdm  4-  ^^wjw^ddm 
+  \UT(bTwJ W Jr  dm  ( i  =  i  ,2,3)  (j  =  1 ,2,3) 

|" \TTWI Wfydtn  4-  J U T/b 7 V  'J IV (f> dsn  1 


\rTWT2Vi)dm  +  \Ur4>Tw'*2W4>dm 
{r*7  W4>dm  4-  J  U  ^V7  iV4.dm 


Mfiq  ~ 

M„n  «  \<l)TWTW<j)dm 


Gn  =  \24>T  WT  Wd)dm 
Kn  =  f <t>r\Vr<t>dm  4-  Jr7cr^x 


-  -  jr7  i^J Wr4>Udm  - 


2  j  r  7 1  vj  IV  (f>  J  dm 


t-  jT7  V/jgdm  +  \gTW44>Uam  ( i  --=  1 ,2,3) 

>;  =  j <j>TWTgdm  -  \4>TWTWfdm  4- 
^  =  1>,1  Fq2  Fl3~\‘ 

The  final  equations  of  motion  are  written  as 


[B  -  28) 

(. B  -  29) 

(3  -  30) 
(£  —  31) 
(B  -  32) 
(B  -  33) 


(3  -  34) 
(B  -  35) 
(/?  -  36) 
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Mqq  «  +  Mqn  U  =  Fq 


Mnq  £  +  Mnn  U  +  Gn  U  +  Kn  U  =  Fn 


(B  -  37) 


{B  -  38) 
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APPENDIX  C  DERIVATION  OF  THE  GENERALIZED  FORCES  OF 
THE  EQUATIONS  OF  THE  MOTION 

To  derive  the  generalized  forces  the  principle  of  virtual  work  is  used.  We  first 
write  the  transformation  between  global  coordinates  and  local  coordinates,  Eqs. 
(C-l)  and  (C-2)  respectively, 

~ i  i  r  i  o  o  iff 

X  =  X0  cos(0)  sin(0)  a:  (C-l) 

Y  Yq  sin(O)  cos(0)  y 

"il  r  i  oo  lr I" 

x  —  ~A"0  cos(0)  —  Yq  sin(0)  cos(d)  sin(0)  X  (C  —  2) 

y  X0  sin(d)  -  F0  cos (6)  -  sin(d)  cos(d)  Y 

From  Eqs.  (C-l)  and  (C-2)  we  can  find  x  and  y. 


w  =  -X0  cos(0)  -  F0  sin(0)  4-  X  cos(9)  4  Ysin(6 )  (C  -  3) 

y  -  XQ  sin(0)  -  F0  sin(0)  4  X  sin(0)  4  Y  cos{fi)  (C  —  4) 

Taking  derivative  of  the  Eq.  (C-3)  and  (C-4)  doing  necessary  substations,  we 
can  find, 

<5x  =  <5  A-  cos(0)  4-  5  Y  sin(0)  (C  -  5) 

dy  —  —bX sin(0)  4  <5  T cos(0)  (C  -  6) 

Using  the  virtual  work  principle, 

5  Wk  =  [  -7'cos(<5  4  O(O))v(O)](<50  4  <5<D(0)) 

4-  [T cos(<5  4  <D(0))]<5jr  4  \T  sin(<5  4  <P(0))](<5y  4  <5v(0))  (C  -  7) 

Based  on  small  angle  assumptions, 
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cos(<5  +  0(0))  =  1 


(C  —  8) 


sin((5  +  0(0))  -  <5  4-  <P(0)  (C  -  9) 

Doing  the  necessary  operations  and  substitutions, 

SU'/(  =  ( 7cos(0)-7<5  sin(0)--7O(O)  sin(0))<5Z  +  (7 sin(0)  +  TS  cos (0)  + 

70(0)  cos(0))<5  Y  +  (  -  7v(O))<5(0)  +  (TS  +  70(0)>5v(0) 


+  ( ~7v(0)<50(0)) 


(C  —  9) 


From  Eq.  (C- 10),  one  can  write  large  motion  generalized  forces  ( Ft ),  i.e.,  Eq. 
(C-I  i)  and  small  motion  generalized  forces  ( Fu ),  i.e.,  equation  (C-12), 


7 cos(0)  -  TS  sin(0)  -  70( 0)  sin(0) 
T sin(0)  +  TS  cos (0)  +70(0)  cos(0) 
-  7v(0) 


(C  ~  11) 


TS  +  70(0) 
-7v(0) 


(C-12) 


APPENDIX  D.  COMPUTER  SIMULATION  CODE 

definitions  of  variables  %%%%%%%^%t%%%%%%%%%%%%%% 

\ 

%  aO-7  -TTiq  coefficients  for  the  rumarical  integration 

%  area  “Cross  sectional  area  of  the  missile  (m**2) 

%  betal  “Nodal  node  ).  % 

%  beta2  ^t-todal  node  2  % 

%  count  -Constant  that  defines  the  element  positions  in  vectors.  % 

%  dtheta  “Desired  pitch  angle (dgrs)  % 

%  dthetadot-Oesired  angular  velocity (rd/s)  % 

j  dt:h«ta2dot=Oesired  angular  acceleration (rd/s**2)  % 

%  delta  “Control  angle (dgrs)  % 

%  E  “Young's  modules (pascal)  % 

%  epsilen  “Parameter  for  integration  coefficient  % 

"a  fn  “Small  motion  force  coefficient  matrix  % 

%  flexifale=Corstant  that  does  the  flexible  part  on  and  off  % 

"s  Inn  -Reformed  small  motion  force  coefficient  matrix  % 

"  fq  “fdtge  motion  force  ’/ector  % 

%  ftima  “Finish  time  of  the  integration  % 

%  force  “External  force (Newton)  % 

%  an  -Gyroscopic  coefficient  matrix  % 

%  grav  -Gravitational  constant (m/s**2>  % 

%  grawee  “Gravity  matrix  % 

%  h"  “Time  step  used  in  the  integration  % 

%  izz  “Horrent  of  inertia  about  the  z-axis(m**4)  ,  % 

%  iota  “Parameter  tor  integration  coefficients  % 

\  kn  -Stifress  matrix  % 

1  k  “fkwfoer  of  the  column  of  the  matrix  % 

%  kp, kv  “The  position  and  velocity  feedback  gains  % 

%  1  “Nunioer  of  the  row  of  tf<e  matrix  % 

1  1 force  “Large  motion  force  vectcr  % 

%  Ini  “Die  length  of  the  missile  (meter)  % 

-  lrm  “Constant  ‘lm*n)  % 

1  mn  -Reformed  small  motion  inertia  and  coupling  coefficient  matrix  % 

%  mgn  “Coupling  inertia  coefficient  matrix  % 

'k  mqq  “Large  motion  inertia  coefficient  matrix  % 

1  ma  “Constant  (mi*area)  % 

l  nphi  “Shape  function  matrix  % 

1  nrphi2  “Second  time  derivative  the  snape  function  % 

i  mu  -Mass  density  of  the  missile (ky/m**3;  % 

%  n  -fkmteer  of  tlie  integration  % 

%  phi  notion  slope  position  I 

%  phidot  -Snail  motion  slope  velocity  % 

%  pnl2dot  -Snail  motion  slope  acceleration  % 

%  rigid  -Constant  that  toes  the  controller  on  and  off  % 

1  rloca!  “local  position  vector  % 
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%  q  “Large  motion  generalized  position  vector  % 

l  cfJot  “Large  motion  generalized  velocity  vector  % 

%  q2dot  “Large  irotion  generalized  acceleration  vector  % 

%  stif  “Stlfness  matrix  % 

%  sforce  “Snail  notion  force  vector  * 

%  slope  “Trajectory  path  slope  % 

%  t«anpfirst=Tenporary  matrix  for  the  Simpson's  rule  % 

%  templast-^Tenporary  matrix  for  the  Slnpson's  rule  % 

%  temp  1 .  • .  % 

%  tempS  *=T<arpcrai.Y  matrices  for  integration  % 

%  tenp7  „ , .  % 

%  tenplO  “Temporary  matrices  for  large  a:id  small  motion  coefficients  % 

%  theta  “Large  motion  angular  position  (pitch  angle-dgrs)  % 

%  thetadot»Large  notion  angular  velocity (rd/s)  I 

%  theta2dot=Large  notion  arvjular  acceleration (rd/s** 2)  % 

%  time  “Sec  % 

%  u  “Small  motion  generalized  position  vector  % 

%  udbt  “Small  motion  generalized  velocity  vector  % 

%  u2dot  “Small  motion  generalized  acceleration  vector  % 

%  un  “Small  motion  position  vector  % 

%  umdot  “Small  motion  velocity  vector  •.  % 

%  v  “Small  notion  deflection  position (m)  % 

l  udot  “Snail  motion  deflection  velocity (m/s)  % 

%  v2dot  “Small  motion  deflection  acceleration (m/s**2)  % 

wrosid  “Residual  acceleration  matrix  % 

%  wdot  'q’ime  derivative  of  the  transformation  matrix  % 

%  wksi  “The  derivative  of  the  transformation  matrix  with  respect  to  the  % 

1  large  motion  generalized  coordinates  % 

%  X  “Large  motion  x-direction  position (m)  % 

%  x  “Local  x  coordinate  % 

%  xc|ot  “Large  motion  x~direction  velocity  (m/s)  j  % 

%  x2dot  “Large  motion  x-dlrection  acceleration (m/s**2)  % 

%  v  “large  motion  y-direction  position (m)  I 

%  ydot  “large  motion  y-direction  ’reiocity (m/s)  % 

%  y2dot  -Large  motion  y-direction  acceleration (m/s* *2)  % 

%  mnnl,knl,  % 

%  kn2,fnl,  % 

%  fn2, fn33,  % 

%  gnl,mqql,  % 

%  fql, fq2  “Tcnporary  matrices  for  the  calculation  of  the  coefficients  % 

%  cl,c2,fl,  % 

%  f2,  f3,  f4“Sliape  function  constants  % 

%  m22,  f22,  I 

%  h2,b2,A,  I 

%  B,F,tanpl4,  % 

%  tenp20  “Terporary  matrices  for  the  controller  % 

%  mll,ml2,  % 

%  m21,m22,  % 

%  fll,hl,bl*=Constants  for  t)>a  controller  % 

%  tolot  “Time  plot  vector  % 
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%  xp.lot  -I.  irga  met ion  x-directicr.  position  plot  vector  % 

%  yplot  “Ij  rge  motion  y-dtrection  position  plot  vector  % 

%  thetaplot“Le  rge  motion  angular  position  plot,  vector  % 

s  xdotplot  -Large  motion  x-direction  velocity  plot  vector  % 

%  ydotplot  “large  motion  y-dl section  velocity  plot  vector  % 

%  thetadctplot  Targe  motion  angular  velocity  plot  vector  % 

%  vplot  =Snvll  motion  deflection  position  plot  vector  % 

%  phiplot  “Sim  il  motion  slope  position  plot  vector  % 

1  vdotplot  “Smell  motion  deflection  velocity  plot  vector  % 

%  phidotplot“Sn  all  motion  slope  velocity  plot  vector  % 

%  Del  tap? ot  -Coitrol  angle  plot  vector  % 

$  dtix»taplot=Oe  ilred  trajectory  angular  position  plot  vector  % 

%  dtlW3tadeg,cteH  adeg,  % 

%  tbetacfcg  “Anc  les  converted  from  radian  to  degrees  % 

%  % 
%%%%%%%%%%%%%%%%%  ¥%%%%%%%%%%%%%%%%%%%%%%%%%^!%%%%%%<%%%%%%%%%%%%%%%%%%%%%%%%%%% 

11%%%%%%%%%%%%%%%',  %%%%%%%%%%%%  MAIN  FROGRAM  *%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%  LEVEL-  1  % 

%  Tills  Is  the  co  (trolling  program  for  flexible  missile.  It  defines  che  % 
-  variables,  dehen lines  the  coefficients  for  explicit  integration  and  mode  % 
%  stap'1  function,  < aj_ls  large  and  small  motion  coefficient  mas  rices  routines, % 
%  large  ard  small  notion  integration  routines,  controller  routines,  output  % 
%  routines  % 

%  % 
%%%%%%%%%%%%%%%*%%%*  <%%%%%ii%%%%%i«%%%%%%%%%%%%%%%%%%%«%%%%%%%%%%%%%%%%%%%%%%%%% 

(X, y, theta, xdot,ydot, thetadot, ln?,betal,beta2,E, izz, force, . . . 
phi, delta, v,  iota, epsilon,  ftine,h, grav,.ra,  vdot,phlclot,n,  Inm,  flexible, . . 
count, mphi,nphi2,sti.f  gravvec,  rlocal,x2dot,y2dot,theta2dot,v2dot,  . .. 
phi2dot]  ™  constants!  j 

(m22,  f22,h2,b2,A,F,  B,  dtl»eta2dot,dthetadot,dtheta,kv,  . . . 
kp,  rigid,  cenpl4, slope,  i:enp20J=corc,tantsla; 

(q,  qiot,q2dot,teip6,wircsid,wdbtt  vksi,  vr, . . . 

1  force,  s force,  um,undot,  fnn,gn,  Ini, mn, count  1, mil, ml 2, m2 1,  f  11,  hi,  ... 
bl ) “const auits2&  (X  y,  theta, xdot,  ydot.  tfvutadot, x2dot,  y2dot,  theta2dot ) ; 

lu,udot,u2dot,  mnnl, knl,Ln2,  fnl,  fn2,  fn33,i»in,gnl,mqn,mqn,  . . . 

t’q,  mqfgl ,  tempi ,  tonp2,  tempi,  renp4,  tenpS,  tesrp7,mqnl,  tenpP.  fq] ,  tenp9,  . . . 

tenplO,  fq2J  -  constcinta2h(v,phi, vdot,phidot, v2cbt,phi.2dot)  ; 

|tplot,xplot,yplot, . . . 

thetaplot, xdotplot,  ydotp  lot,  tMtadotplot, vplot, j>hl plot,  vdot  plot,  . . . 
phi  riot  plot,  the  .aritea.  delbyplot,  •  lthetaplot,dthetacfcy,  deltadeg J  =. , . 
conshants3  (fti.no,  h) ; 

(aO,  al,  a2,  a3,  a 4,  a5,  a6,  a  / 1  »•  ccoof  iepsilon,  iota,  h) ; 
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Icl,c2,  fl,f2,f3,  t4)”«phAConsUm/betal?beta2) ; 
for  time  •*  O.Oshiftiroe, 

[w, wdot, wksi, wresid. lforce, sforce,  urn,  umdot)  —  wmatderiwai  iant  i , . » 
theta,  X,y,xdot,ydot,  thetadot,  forex*,  phi,  ctelta,v,vdot,phidot,... 

wresid,wdot,wksi,w,  1 force,  s  force,  urn,  undot} ; 

[mqq,nryi,  fq,  fq2]—lrgcof  (wksi,  rlocal,rrphi,ixti,  w,  1  force,  wresid, . . . 
wdot,  uindoc, q rawer,  ma.  lm,  n,  mqq,  mqn,  fq,rtK^l,tenp7,mqnl,  . . . 
terrp8,  fql,  tu.p9,  tenplO,  lwi,  cl, c2,  f  1,  f2,  f3,  f4,betal, beta2,  fq2)  ; 
if  flexible— 1, 

(rn,gn,  kn,  fnnl-amlcof  faphi,  w,mqn,mqq,  fq,wdot,wresid,irphi2,  stif,  . . . 
grawec,  sforce,ma,  lm,  n,  lmn,  fnn,gn,kjn,mn,nnnl,knl,kn2,  fnl,  fn2,  fn33,  . . . 
wui,  rlocal,d,c2,  fl,  f2,  f3,  f4,betal,beta2)  ; 

(u,udot,  u2dot)=*intsml  (aO, al, a2,a3, a4,a5,a6,  a7,u,udot,u2dot, . . . 
t«Tpl/tenp2,t»3Tp3,  terp4,tenp5,nii,gn,  fnn,  kn) ; 

end 


[q,qdot,q2dot]  -  intlrg  (h,  aO,  a3,  a6,  a7,  mqn, mqq,  fq,  q,  qdot,q2dot, . . . 
tatp6,  u2dotf ; 

(X, y, theta, xdot,ydot, thetadot, v,pbi, vdot,phidot,  x2dot,  . . . 
y2dot,  tlieta2dot,  v2dot, phl2dot ] -chvar (q, qdot, q2dot,  u,  udot,  u2dot ) ; 

if  rigld=l.- 

I delta, dtheta, t«rp20]-rigidcontrol  (mcjq,  fq2,  force,  theta,  mil, ml2,m?l,m22, 
fll,f22,hl,h2,bl,b2,A,!r,B,dtheta2dot,U)etacbt,dt)ietadot,  • .  • 
dtheta,  kv,  kp,  tenpl4,  time,  slope,  tenp20) ; 
end 

(tplotfjqolotjyplot,  thetaplot,  xdotplot,  count,  deltaplot,  dtlmtaplot)  -gplot  ( 
court,  tiire,X,  y,  tlveta,xdot,  tplot ,  xpiot,  yplot,  thetaplot,  xdotplot,  thetacteg,  . 
delta,  deltaplot,  dtheta,  dtlietadeg,  dtliet aplot ,  deltadeg) ; 

(ydotplot,thetadotplot,  vplod,phlplot,  vdotplot,  phldotplot,  count  J,=gplot2  < . . 
count,  ydot,  thetadot,  v,phi,  vdot,phidot,ydotplot,  tl>etadatplot,  vplot,phiplot 
vdotplot , phidctplot ) t 

count 1-count 1+ 1 

if  count.  1—0,  . 

fcroe-0; 

end 

if  countl— 0, 

tl«X; 

t2=y; 
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t3=theta; 

t4=xdot; 

t5~ydot; 

t6=fctaetacbt; 

t7=v; 

t3=phi; 

t£=vdot; 

tlO-phidot; 

tll=x2dot; 

tl2«=y2cfot; 

tl3=thet.a2dot; 

tl4”v2dbt; 

tl5=phi2dot? 

tl6=«delta; 

clear  X  y  theta  xcbt  ydot  thetadot  v  phi  vdot  phidot  delta 
clear  fnl  fn2  fn33  fnn  fq  fql  gn  gnl  kn  knl  kn2  1 force  nr>  imn 
clear  nnnl  irqn  mqnl  mqq  mgql  q  qdot  q2dct  sforce  tenpl  te*ip2 
clear  te«p3  tenp4  tenpS  tetp6  tenp7  terrp8  tetp9  u  udot  u2dot 
clear  im  undbt  w  wdot  wksl  wresid  count  1  tenplO 
clear  x2dbt  y2dot  theta2dot  v2dot  phi2dot 

clear  mil  ml2  m21  fl  hi  bl  teirpll  terpl2  tenpl3  tenpflrst  terrplast 


X=tl; 

y=t2; 

theta=t3; 

’  xdot=t4; 
ydot*=t5; 
tiietadot“t6; 
v=-t7; 
phl«t8; 
vdot=t9; 
p!ildot=tlO; 
x2dot=tll; 
iy2dot=t.l2; 
Uieta2dot**tl3; 
('2dot=tl4; 
ptti2dot*=tl5; 
ctelta«tlo; 


|q,cjdot,q2dct,tenp6,  wresid, v*jbt,wkai,w, . . . 

1  force,  sforce,um,  undot,  fon,gn,kn,nn,countl,mll,ml2,m21,  fl,hl,  . . . 
blJ-con3tants2&(X,y,  theta,  xdot,ydot,thetadot,x2dot,y2dot,  . . . 
theta2cbt) ; 

lu,  udot, u2cbt, nnnl,  knl,  kn2,  fnl,  fn2,  fn33,  iron, gnl,  mqq, mqn,  . . . 
fq,mqql,  tenpl,  tenp2,  t.errp3,  tcnp4,  tenpS,  tertp7,mry«l ,  teip8,  fq!,tenp9,  . . , 
tenplO]  -  oonstents2t>  (v, phi,  vdot, phidot,  v2dot,phl2dot.)  ; 
end 
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lplotl,plot2, plot3, plot4, plots, plot 6]  =allplotl  (tplot,xplot,yplot, . . . 
tlvetaplot,  xdotplot,  ydotplot,  tlietadotplot,  dthetaplot) 
tplot7,plot8,plot9,plotl0,plotll]=allplot2  (tplot,  vplot,phiplot, . . . 
vdotplot,  phidotplot,  cfeltaplot) 


.  LEVEL  2  , 

l  -ruis  function  determines  constant  values  and  initial  values  of 
%  variables  and  matrices. 

% 


f iction  (X,  y,  theta,  xdot,  ydct ,  thetadot ,  lexSie,'  *. , 

ptii2dot]  ■*  constants  1 


area=0.0113097; 

1111=4; 

be'  al-=4 . 71 2398/lm; 
beta2=7.853981/lm; 
count=0; 
delta=0; 

E--’2.0e+ll? 

e{isllon=0.25; 

f lexlble=l; 

forct=30000; 

ftin»e=0.4; 

grav=-9.9066; 

li=0.001; 

iota=0 .5; 

izz=0. 000010178; 

HM-7861..05; 

tv-10; 

Jn«i=Wn; 

tin  area*nu; 

H>i=o; 

fliidbt*0; 
t  tir*t  a  pl/4; 

t  lietadot-O; 

X  0; 
x<  lot-0; 

\r=0; 
wkvt-O; 
ydot“0; 
v  <->; 

ti^hi*  zeros  (3,2) ; 
i;|'hi2=zeros  (3, 2) ; 
;,t  if-zerOs  (J) ! 
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a»  <*»  & 


stif  n,3)*45*izs? 

giaw^-zeros  (3, 1) ; 

grawac  (3)*egrav; 

rlocal=»2eros  (3,  i) ; 

x2dot*0; 

y2dot»0; 

theta2dot“0; 

v2dot-0; 

phi2dbt-0; 

%  % 

%  This  function  determines  constant  values  and  initial  values  of  the  * 

*  variables  and  matrices.  % 

%  % 


function Jm22,  f22,  h2,b2.  A,  F(  B,dtheta2dot ,  dthetadot,  dtheta,  kv, . . . 
kp,  rigid,  te«npi4,  slope,  temp20]  “Crxistantsla 

itv'2-d); 

f22=0; 

b2«0; 

VO; 

F  0; 

BO; 

dthet«2dot"0f 

dthetacfctm); 

dtlK»ta«0.7854;%45  Degree 
kv=6; 
kfr  9? 
i  jj  g  Id*-  1 ; 

to(pl4«zeros  (2) ; 
s .  ope**"  2  r  .  *'>• 
teaip20=0; 

vimmummummmmmmmmmMwmunonumnnmmu 
1  % 

'  T»ila  function  determines  constant  values  and  Initial  values  of  the  % 

%  variables  and  matrices.  I 

*  % 

imntmmiutktttmmmummimmuiVdtmmmimmnmuui 


function  (cj, cpfcit,  g2dot,ten|)6,wresld,wdr>t,wksl,w, _ 

1  forc-e,sforce,um/umdot.,  f im, gn, lui.mn, ctXMitl.mil, ml 2, m2 1,  fll.hl,  .  . . 
bl  J  -constants2a  (X,  y,  theta,  xdbt ,  yrt* ,  thetadot,  x2dot ,  y2dot,  theta2ctat> 
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wm 


WWMWP aPWM 


q --zeros  (3, 1); 
qdot- zeros  <3, 1) ; 
q2dot»zerce <  3, 1 ) ; 
teip6-zercs (3, 1) ; 
wresich-zeros  (3)  t 
wdot*=zeros  (3) ; 
wksi-zettw  <3, 9) } 
w=zeros  (3) ; 
lforce-zeros (3, 1) ; 
sforte-zeros <2, 1) ; 
unp>zeros(2,l)  / 
umdot-zeros  (2, 1)  / 
fnn-zeros (2, 1) ; 
gn^zeros (2) ; 
kn -zeros (2); 
mn-zeros(2)i 
count  1-0; 
q  ( 1)  >=0C; 
q(2)~y; 
q(  3)  -theta; 
qrlot  (D-xdot; 
qdot(2)«ydot; 
gptot  (3) -thetadot; 
q2dot(l)-x2cfc>t;  r 
q2dot (2)-y2dot; 
q2dot  (3)-theta2doc; 
mll=zeros(2) ; 
ml2=zeros  <2, 1) ; 
m21=zero3  (1, 2) ; 
f 11 -zeros (2, 1); 
hl=-zeros(2,l); 
bl- zeros  (2, 1) ; 

intiHttmmmmwmmiuimmwmummHummtuimnitiu 


%  % 

*  This  function  determines  constant  values  and  Initial  values  of  the  % 

%  variables  and  matrices.  % 

%  1 


immimumiuumHmimimmmmmuutmmmmmiumm 


funct  ion  fu,udot,u2dbt,  avail,  knl,kn2,  fnl,  fn2,  fn33,s*v»,gnl,  wqq,»aqp,. _ 

fq.  mcf^l ,  twpl,  temp2,  t«arp4,  terrp5,  t«Rp7,mqnl,  ten^e,  fql,  tenp9,  . . . 

toiplO,  fq2 )  »  constant  a2b(v,pfil,  vdnt.phidot,  v2<**.,phi2dDt) 

u-zerc*s(2s  1) ; 
udbt -zeros (2,1); 
u2rtot— zerna  <2, 1)  ; 
mini -zeros  (2) ; 
kill -zeros  (2) ; 
kn2-zeros (2) ; 
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fnl-neroa 12, I) ; 
fn2«zeros(2,l) ; 
fn’3“zerca (2, 1) ; 

~n>n -zeros  (2)  ; 
c,nl-«zeifXKj(2)/ 
nr-irzeroa  (3) ; 
in<?i~Eero»  (3,  2)  y 
fq=*zeroa(3„  11/ 

mqql«**eros(3)/ 
tenpl»zeixe  <2, 1 } ; 
tersp2«zeras  (2, 1) ; 
teropB-zeras  (2,1); 
teipfl.  ieros  (2, 1)  ; 
t«tpS~zeros  12, 1J ; 
tesFp7»zerm  O) ; 
mqnl-zerca (3, 2) ; 
t^rpfl-zeros  <3)  jr 
tql-«»r©<*<3,l)# 
fq2»zeios  (3, 1) s 
tatp<>-zeto3  (3?  ? 
fc  «?:>?.  10=zr*roo  (3)  ; 

u(li«v; 

u(2)«phi; 

ucbt(l{-vdotf  r 

ucbh  (2)»phidbt; 
u2dot  ( 1 }  «v2tfc}-  ? 
u2dot  (2)  ~yhi2rtbt7 


%  % 

%  Tills  function  detercdnea  the  Initial  values  of  the  outp-it  matrices.  % 

%  \ 


function  (tplofc,>qpiot,^>Iot, — 

tlietaplot,  wrfc*f>lot,  yciotplot,  thetacJotplot,  vp  lot,  phi  plot,  vrlotpiot,  . , . 
phktotplot,thet®deg,tfait^lof»<*he*:«plo*,ci±*tacSeg,ciBitadeg)  «... 
cofMt.«nts3  (f  tlfRO,  h) 


tpl ot-xeros 1 1 ,  (ftiw*s/h)  tl)  j 
xpfot-zerosU,  (ffc.iw./h!  +  1) ; 
yplot-zerce  i  1,  (ftin«/hi  *1) ; 
tlietaplofc-Twrow 4) ,  (ftiow/h)  *1)  s 
xclotplot“iMiro8  ( i ,  i ft  laie/h)  ♦  1 1 ; 
ydotplc*. -zeros  (I ,  tf  tims/)i)  +11  j 
tlietackitplot -zeros  U,  (Cthme/h)  +  i ) ; 
vpl ot  -zeros  1 1,  ( ft iec/h)  » 1 ) ; 

phlpioc-xeroo;!,  (ftimei/h)  +11? 


vdotf>lot-r.6EOS  \ ,  'ftime/hS  +1)  ; 
phidotplot-zerot  1,  (ftime/h)  +1)  ; 
tlietadegK); 

daltaplot-zerc:  tine/h)  +1) ; 

dthetaplot-zerc  ftijne/h)  +3.) ; 

dtlietadeg-O; 
daltadeg-Oj 

%  % 

%  Formulating  the  coefficients  used  in  the  sequential  integration  method.  % 

%  % 


function  (30,31,32,83,84,35,36,371  -  ccoef (epsl Ion, iota, h) 

20*1/ (epsilon* (h*h) ) ; 
al=iota/ (epsilon*n) ; 
a2*l/ (epsllon*h) ; 
a  3---0 . 5/epsi  lon-1 ; 
a4=iota/epsilon-l; 
a5-h/2*(a4-l); 
a6*h*  (1-iota) ; 
a7=iota*h; 

* 

mmnmimimmuvmnnnmutnmmimmmwniiHunimm 

%  % 

%  Formulating  the  coefficients  used  in  the  mode  shape  function  matrix.  % 

*  % 


function  (cl,c2, fl, f2,  f3,  f4)mnptiicons (1m,  petal,  beta2) 

cl  >  (sin  (beta  1* Ins)  -alnh  (beta! *lm) )  /  (-cos  (betal  *Lm)  4 cosh (betal»lm) )  ; 

c2-  (sin  (bata2*lm)  -alnh  (beta2‘lm) )  /  (-cos  (bet«2*lm)  rcosh  (beta2*lm> ) , 

e*2*beta2-2* (c2/cl) »betal; 

f 1-1/ (2*cl) * (c2*betal) / ( (cl'2) *e) ; 

f  2—  (c2/  (cl*e) ) ; 

f3~-  (betal/  (cl*e) ) ; 

fi-l/e; 


i  % 

'*  This  fiEKtlon  calculates  the  matrices  which  changes  with  time.  % 

*  % 


function  (w,  wdot,  v.ksi,  wresid,  i force, sforce,  tm, undot  j  . . 
wmatderi  vw.r  iant  (theta,  X,  y,  xdot ,  ydot,  thet adot,  force,  pi  li,  del  ta,  v, .  „ . 

T\:L't ,  Ubfc ,  wresld,  wdol'  ,  wksi.,,  w,  1  force,  3  force,  un,  uirJot ) 

%  Conp'Utakion  of  fctje  transformation  iratrix  w  % 

m%mt%%m%%%%%%%%%n%%%%u%%%%%%%%%%%%%%%%«%mm%%%%%t%%%%««%%%%%%m%%<k 

w(lf  ,l)«l; 
w(2,l)=X; 

w  (2, 2)  “coo  (theta; ; 
w  (2, 3)  **-'•'  (theta) ; 
w(3,  l)-=y. 

w  (3, 2)  ’=sin  (theta) ; 
w  (3,  J)  -cos  (theta) ; 

«%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%  Computation  of  the  time  derivative  of  the  transformation  matrix  WDOT  % 
%%%<&«%!*%%%%&%%*%%%%%%%%%%%%%*%%%%%%%%%%*%%%%«%%%%%%%%%%*%%%%%%%%%%%%%%%%%%%% 
wdet (2, l)=xct>t; 

wdot  (2, 2)  '=-thetadot*sin  (theta) ; 
wdot  (2, 3)  ~--thetadot*cos  {tiieta)  j 
wdot  (3,  IKydot; 

wdot  (3, 2 )  ^thetadot  *cos  (t)ieta) ; 

wdot'0, 3)  — thetack>t*sin  (tiieta) ; 

/ 

%  Cornputation  of  the  derivative  of  the  transformation  matrix  with  respect  % 
%  to  X,  Y,  theta i  WKSI  % 

wksi  (2, 1)**1; 

wk.ji  (3,  4)  =  *.; 

wksi  (2,  h)  -»-sln  (tiieta) ; 

wksi  (3, 0)  -cos  (f lieta) , 

wksi  (2,  9)  —“-cos  (tiieta) ; 

wksi  (3,  9)  “-sin  (tiieta) ; 


»%%%%%%%%%%%%% 

%  Cdnpfitation  of  the  the  resickia)  acceleration  nvitrix  WKES1D  % 

ummtvmumuuiimtniHiuutunmuHununuHtnimmmm 

wtesld(2,2)  —  (thetadot~2>  *cu3  (tiieta)  ; 
wresid(2,  3)-(thetacbt'v2)  *sln  (tiieta) ; 
wresld(3, 2).—  (thetadot^)  *sin  (theta) ; 
wresid (.7, 3) •—  (thetactotA2)  *cos  (tiieta) ; 

imiuvuiiHiiuHHHHimtniHimmmmHunmutimmimmui 

1  Computation  of  ttie  larxje  motion  force  nvitrix  % 

'tnuiuniunitimmniimuuitttmtmuitnmmmtmt’iintniim 

]  force (1) -force*  (cos  (tiieta)  -delta*sin (theta)  -phi*sin  (tiieta:  )  ; 
i  force  (2)  -force*  (sin  (tiieta)  tdeJtta*oos  (l)ieLa)  iphl  *cos  (tiieta) )  ; 
ifocoe  i3)“~ioi:oa,*v; 
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%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%  Confutation  of  tlie  small  motion  force  natrlx  % 

%%%%%%%%%%%%%%%%%%*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
sforce  (1}  -force*  (phi+delta) ; 
sforce  (2)  **v*forca; 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%  Computation  of  the  snail  motion  position  vector  % 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

um(li'=v; 

un(2)«phi; 

%%%%%%%%%%%%%%*%%%%*%%%%%%%%%%%*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%  Computation  of  the  small  motion  velocity  vector  % 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%£%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
vundot  (l)=vdot; 
uittiot  (2)=phicbt; 


%  This  function  calculates  values  for  the  large  motion  coefficient  matrices% 


function  lmqq,rnrpi,  fq,  fq2]“lrgcof  (wksi,rlocal,nphi,i«n,w,  lforce,  . . . 
wresid,  wdot,  unrfot,grawec,ma,  lm,  n,mqq,mgn,  fq,inqgl,  tenp7,mqnl,  . . . 
tonp8,  fql ,  temp9,  tenplO,  lmn,cl,c2,  fl,  f2,  f3,  C4,betal,beta2,  fq2) 

%  Creating  ttie  coefficient  matrix  for  MJQ  % 

for  i=l:3, 

terTpl-wksi{:,3*(i-l)tl:3*i); 
for  j=l:3, 

teip2=vksi(:,3*<)~l>  O  s3*  j) ; 

taip7=>tenpl 1  *tenp2; 

rnqql  (i,  j) -“quadrat  ('astringll',  1,1, n,  Imn,ma,cl,c2,tetal,beta2,  fl,  f2,  f 3,  . . . 
f  4,nphi,0,  0,  clocal.  Carpi,  t3:p2.  *sn,  0,0, 0,  0,tenp7, 0,  0,  . . . 

6,0,0); 

raid 

oik) 

nr  f  [  TiVi*  (mqqi)  ; 

%  Creating  tire,  coefficient  matrix  for  HJi  1 

numumMmHmmnmuuummummmmnmimimmnnn 

for  1-1:3, 

t  «rp>^»ksi  ( : , 3* (i-1)  *1:3*1)  ? 
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teipO^erp '  *w; 

mqnl  (i,  s  )“cjuadkTiat  ( * astring22 ' ,  1, 2,  n,  lmn,ma,  cl, c2, betal,beta2,  fl, £2, £3, . . . 
f  4,  nphi,  w,  0,  rlocal ,  0, 0,  urn,  t«rp,  0, 0, 0, 0,  texp8, 0, 0, 0, 0) ; 


end 

mqpFH0*  (rttfil) ,* 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%'*%%%%%%%%%%%'*%% 
%  Creating  the  coefficient  mayrix  for  FQ  % 


%Caiputation  of  the  FQ  coefficient 
for  i=i:3, 

tarf>=wksi  { :,3*  (i-1)  +l:3*i) ; 
tenpSM:eiTp'  *wresid? 

f  Jui)^Sknat  Castring33' ,  1,  l,n,  lWn,ma,  cl,c2,b=tal,beta2,  f  1,  f2,  f3,  f4, 
r phi,  0,  wtesid,  rlocal,  0, 0,  um,  tenp,  wdot,  undot,  gravvec,  0, 0, tenp9,  . .  - 
tenpl0,0,0)  ,* 


aid 

fq2=ina*fql; 
fq-fq2+-l  force; 


%  % 

%  This  function' calculates  values  for  small  notion  coefficient  matrices.  " 

%  % 


function  (mn,gn,kn,  fnn)=-smlcof  (nphi, w, rtcji, mqq,  fq,v*Jot,wresid,nphi2,  •  .  . 
stif,grawec, sforce,nva,  lm,n,  lmn,  fiui,gn,  kn,nri,iminl,knl,  — 
kn2,  fnl,  fn2,  fn33,tram,  rlocal,  cl,c2,  fl,  £2,  f3,  f4,Lx?tal  ,beta2) 

%  Creating  ttie  coefficient  matrix  Ml 

rwui=quadnat  (’astring44’,2,  2,n,  In*i,ma,cl,c2,betal,t>eta2,  f  1,  f2,  IJ,  f  4,  ... 

n|3ii ,  w,  9, 0, 0,  C,  0, 0, 0, 0, 0,  0, 0, 0, 0, 0, 0)  ; 
n»u.l-ina*n«ui; 

irmmuil-tngn'  *  Inv  (mqq)  *n£jn; 

%f %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%« 

%  Creating  tlie  coefficient  matrix  CN  l 

inmuiKiimmummmnutmumimummnmnmnuHmnn 

gnl*sguactnat  1  'astring55' ,  2, 2,  n,  ln«i,ma, cl, c2,l3etal,beta2,  f  1 ,  f 2,  f 3,  f  4,  ... 

mill,  w,  0,0, 0,0,0, 0,  wdot,  0, 0, 0, 0, 0,  0, 0, 0) ; 
gn=ma*gnl; 

nuumHmvmimmmmnmmummmmmnmimnmmimi 

%  Creating  t!ie  coefficient  matrix  KN  1 

knl*=quadmat  ( *a3tring66',  2, 2,  n,  ln»i,ma, cl,  c2, fetal , beta2,  f  1,  f 2,  f  3,  f  1,  ... 


80 


nphl ,  w,  wre3id,  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0) ; 
kn2-quadknat  Castring77 ’ ,  2, 2, n,  luti,ma,  cl,c2,betal,beta2,  £1, 12,  f  3,  £4,  . . . 

0, 0, 0, 0, 0, 0, 0,  0, 0, 0, 0, 0, 0, 0, 0,  nphl 2,  stif ) ; 
ki  *=ma  *kr;l+kjv2; 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%'&%%%%%%%%%%%%%%%%%':. 
%  Creating  trie  coeCficient  matrix  EN  " 

fnl-quacknat  Castri.ng88', 2, 1 , n,  In*i,ma,cl,c2,betal,beta2, fl, £2,  f3,  f*t,rcphi,wf . . . 

wresid,  r local,  0,  0, 0, 0, 0, 0,  grawec,  0, 0, 0, 0, 0, 0) ; 
tn33=ma* (fnl)  +  sforce; 
fnn=fn33-nqn’*inv(mqq|)  *fq; 

%  % 

%  This  function  iirplicitly  integrates  the  small  motion  equation.  % 

%  » 


function  (u,  U'Jot, u2dot)“>int.,>inl  (a0,al,a2,a3,a1fa5,a6,  a?,u,udot,  . . . 
u2dot,  taipl,  tatp2,  tenp3,  tatp4,  teep5,wi,gn,  £nn»  kn) 

taipl“s0*uia2$pdot>a3*u2dot? 

tenp2>=a  l*u+a1*  ticbt  f  a5  *  u2dbt ; 

te»p3=ann*taTpl; 

terp'I^Ti*  tatp2  ; 

tenpl-f  nn+tenp3+ taipl  ; 

t«ip2=u; 

tenp5»kn+a0*mi+al*gn; 
t  en  p  1  «*t  enp5  \  t  en  p  1 ; 
u=taipl; 
teip.3=u2dot ; 

u2dot-a0*  (u-tenp2)  -a2*udct :-a3l't.enp3; 
udotaxiat  *a6‘ terpJ* a7 * u2det ; 

% 

%  Ibis  function  explicitly  Integrates  tl>e  large  motion  equation. 

% 


function  [q,  <jiiot,qf2(lot  J”lntlrg  (h,aO,  a i,  36,  a  1  ,t»  ,n,ntf  ],  fq,  q,  qkil ,  .  .  . 
f|2dbt ,  te  ip6,  u2dot ) 

q  r|*h*qdot+a3*c}2ctot/a0; 

<  f  kit  =« plot  +a6  *q2ck 
t  eii|>6“t6tiy>6  *  mt  jn  *  u2ri>t  i 
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fqn=fq~fc€a<pb; 
q^dot^lnv  (inqcj)  *fqfi; 
cy.iot“qdot  4q2i±>t  *a7  ? 
q^q  Kj2cbt/a0; 


% 

%  Tills  function  updates  tl>e  c.itre  dependent:  values. 

% 


f'irsctloii  5X,y,  theta,xdr.t,y'dot, fclietadofc,v,plil,  vc>jfc,pldrJtotf  x2dot, . . . 
y2dot,  theta2dbt,  vAtot,  phi2dofc.]°chvar  (q,  q-Jot,  q2dot,  u»  udofc,  u2dot) 

X=q(.I); 

y--q(2); 

Mieta^qO) ; 
xdbt«cy:iot  (1)  ? 
ydot^qdbt  (2) ; 
tJ>etadot--qpot  (3)  / 
x2dot»q2ciot  (1)  ; 
y2dot=q2ctot  (2) 
tlieta2dot*xj2dot  (3) ; 
v“=u(l.>; 
pM=u(2); 
vdot“udot  (1)  ; 
phldot-“udat  (2)  ; 
v2cbt«u2dot  (1* ; 
phi2dot=u2dot  (2) ; 


4  Ill's  function  controls  rigid  or  flexible  missile  with  riqld-body 
%  coi  it  roller.  lire  control  angle  (delta)  is  limited  to  10  degrees. 

% 


function  (delta,  dtlreta,  t@ip20)-rigictecmtrol  !ntf f  <\7,  force,  l  beta,  ml  l,nl2,  .  .  . 
m?l,m22,  fll,  f 22, hi,  h2, bl , b2.  A,  if, B, dtlieta2dot,  tlietadot,  dtlietadot ,  . . . 
dtlreta,  kv,  kp,  tenpll,  time,  slcj  <e,  tenp20> 

if  (time  >»  0)  &  (t.lit»  <  0.1J, 
dtl»eta2<A>t=0; 
dtlietadot-O; 
dt(  retard/ 4; 
tenp20-dthet  a ; 
elseif  time  >»  0.1, 
idti»eta2db(;»0; 
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dtl  tetadofc=s  lope; 
dtlieta~tenp20+3lqpe*  (time-O.l) ; 
an 

mil  (1,1)  ■mg?  (1,1)  / 
mil  (1,2)  "tot? (1,2) ; 
mll(2,  l)“mqn(2,l) ; 
mil  (2, 2)*urjq(2,2) ; 
ml2(l)=mqq(l#3); 
ml2(2)=mfj(2,3); 
m21(l)*=nKR(3,l); 
m21  (2)  men  (3,2*; 
nt22=mqq(3,3) ; 


£ll(l)-fq2(l); 
fll(2)-fq2<2); 
f22“£q2  (3) ; 

hi  (l)”forcc*co3  (theta)  ; 
hi  (2)  =force*siir«  (theta) ; 

bl  (l)«~focce*sln  (theta)  ; 
bl (21 =fotce*cos (thets) ; 
taipl4=inv(mll)  ; 

Atr\22  -m2 1  *  teipl  4  *ml2  ; 
e=f22-m21*teipl4*  (fll+hlf ; 
B=an21*tenpl4*bl; 


delta=-inv  (B)  *  (A*  (dtheta2cL>t-kv*  (thetadot  -dthetadot )  -kp*  (tlieta-. . . 
dtlieta) )  -E ) ; 

if  delta  >=■  0.174532  ,  %10  degree 
delta^O. 174532; 
end 

if  delta  <-  -0.174532, 
delta— 0.17  4  532; 

elself  (delta  <  0.174532)  I  (delta  >  -0. *74532), 
delta-delta; 
e«*3 

nunuimmimimumummtiimmmnuimnvmnmmmitm 

%  \ 

*  71. is  fufvction  calculates  the  elements  of  matrices  lot  oulfxit  % 

%  % 

mmntmuimimmmimmiimummnmmnmmmuvimm 


function  |tplotr  xplot,  yplot,  thetaplot,  xdot plot,  count  .delt  aplnt ,  . 
dtlietaplot]  “gpiot  (count,  t  ime,  X,  y,  tlieta,  xdot,  tplot .  xplct .  yplot,  .  . 
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tbetaplct,Mdbtplot,  tbstadlag,  delta, ck^Lfcapiotfdthetajdthetacleg, . . . 
dtlietaplot,  daltacieg) 

count.«ccMnt+ 1 
tSietadeg*  '180*thet&.])  /pi; 
dthetacteg™(180*dtheta) /pi? 
deltadey- (lRO*delta) /pi; 
tplot  (count)  *=tlme; 

Mplot (count) *X; 
ypiot (count) «=y; 
ttietaplot  (cormtS^Thetacteg; 
xdotplot  (count)  **xciot( 
daltaplot  (count)  ^.leltacjeg; 
dtl>dtap.lot  (count )  “dthetadeg; 


%  % 

%  Tnis  function  calculates  the  elements  of  matrices  for  output  % 

%  % 


function  (ydbtplot, tlietadotplot, vplot, phi plot, "/dot plot , ph  Idol: pi ot , count  1  «■ — 
gplotZ  (count,  ydot,  fchetadot,  v,phL,  vdot,phictot,yctot.p.lot,  ttietucv.^plcL,  .  . . 
vplot,  phipiot,  vdotplot,  phldotplot ) 

ydotplot:  (count)  =ydot; 
tliet.adotplot  (count.) “thetadot; 
vplot  (count)  “v; 
phiplot  (count >  “pl'ii ; 
vt  lotp)  ot  (cthji  it )  'vk't  ; 
phldotplot,  (count )  ‘pliAdot ; 


luumumttwmnHnmtHmthmmHnusHKmmummuiH'dt 

%  * 

4  This  function  plots  tlie  output  graphs.  % 

%  % 


function  (plot  1,  pi ot2, plot 3, plot 4, plot5, plot 6| -a  1  lplotl  (tplot,  sqplet,  yplot,  . 
lift  apl  >t ,  xtA.it  pi  ot.,  yvint.pl  ot ,  i.tiet.adbiplot ,  dthetaplot ) 

plot  1  plot  (tplot, iqilot  ) , pause 
plot  i  p  ot  (tplot  ,  yplot) , pause 

pi  ot  1-plot.  (tp.lot,  tbetaprU ,  •  tplot,  dthetap*  at, ’ t  ’) ,  pause 
plot  4-plot  (t  .plot,  mint  plot) ,  pause 
plot  5-plot  (tplot,  ydotplot: ) ,  pcuse 
plot  '.“plot  (tplot.,  thetadot  pi  :>t ) ,  pause 


mn 


ua 


%%%%%*&%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

%  % 

%  This  function  plots  the  output  graphs.  % 

%  ~  % 


funct ion  (plot7, plots, plot9, plotlO, plot 11)  >=al lplot2 (tplot , vplot, phiplot, _ 

vdotplot,phidotplotf  del taplot ) 

plot7^  plot (tplot, vplot) ,  pause 
plotS^plot (tplot,  phlplot) , pause 
plot.9plot  (tplot,  vdotplot) ,  pause 
plotlO=p>lot  (tplot,  phidotplot) ,  pause 
plotll-plot (tplot,  deltaplot) , pause 


%  level  3  % 

%  This  function  will  integrate  the  coefficient  matrices  of  small  and  large  % 
%  motion  over  the  length  of  the  missile  using  Simpson's  integration  metiiod.  % 
%  % 


function  aa  •»  guadint  (F,  1,  k,  n,  int,  jk,  cl,c2,betal,beta2, . . . 

Cl,  f2,  l  3,  f ‘i , m.hl , w,  wres i d,  x~ local, teny>l,temp2,um,  tenp,wdot,»mdot,  . . . 
grawtx;,  t.amp7,  tenpB,  t«*ip9»  tenpl0,nptii2, 3tlf ) 


tempi  1  zeros  (l,k*n) ; 
twpl',"7.ero8(l,k) ; 
t  erp  1 3  ■-  z eros  f  1 ,  k ) ; 
t  arnpf  i  i  ist.-y.eioe  ( 1 , k) ; 
templiv,! .«-ieios(l,k) ; 
«a*zerr«  (1,  k) ; 


k-0; 

far  i-  1  :n, 

X“X 'int; 

itirpll  (l-l)  +  l:  l*k)-f«rvU  (F,x,cl,c2,betal,bet«?,  fl,  f2,  f3,  t4,  . . . 
spin ,  w,  wi.es ' d,r  local,  tempi,  t  mp2, um,  tenp,  r local ,  wdbt.um Sot.,  . . . 

<  j  i  a  v  vec ,  t  enp7 ,  t  ttspB ,  t  enp9 ,  t  &  spl  0,  n  j  hi  2 ,  at  i  f )  • 


a  d 

for  rrr~  1  :k, 

for  ,  *»  2:2:  (n-l> , 

oin)l2{:,m)  «  tfsrpl 2  ( : , m)  ♦  t  enpl  1  ( : ,  k*  (i-l)tin)  ; 
«i 
ail 
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for  p---  l:k, 

i  or  j  =  3:2;  fn-1) ,, 

te«pl3(:,p)  -  tajpl3(;,p)  +  terpll':,  K* (  j-1)  +p)  ? 

aid 

end 

tr  ipf  1  rst  ( : ,  1 :  k)  --tempi  1  ( : ,  1 :  k)  / 
terplast  l:k)=tarpll  (:,n*k-(k-l)  :r.*k) ; 

fi.T-  (Jrit/3)  *  (tenpfirst+4Atenpl3i2*tejrpl2+tefTplast)  ; 

%  % 
%  IT. *3  function  creates  tiie  product  of  matrices  for  integration  of  M3Q 

%  coefficient,  matrix.  % 

% 

<&%%%%%%%%«%%%%%%%%%%%%%%%*,%%%%%% 


function  yll^astringll  (x,  cl,  c2,  Petal,  beta2, €1,  f:2,  f3,  f4,nphi,w,  . .  . 
wi es id,  rlocal,  tempi ,  tenp2, um,  temp, r  local,  wdot ,  undot, gtawec,  tenp7, . . . 
t .  upfl,  tarp9,  tempi  0,qj  hi2,  st  '  t) 

117  hi  (3,l)“fl*  (cl*  (cos  (beta*  *x)  icot.ii  (to* :a 1  *x) )  fsln (beta)  *x)i . . . 
s  i  nh  (beral  *x‘ )  i  f  3*  (c2*  (cos  (bela2*xii  teosh  (beta2*>.) )»•... 
i n  (bet.a2*x)  isinli (beta2*xl )  ; 

nv lit  (3,  ?.)  -f2*  (cl  *  (coa  (total  *x)  teosh  (total *x) )  i  sin  (net  a  1 *x) t . . . 
r.inh  !total*x) )  til*  (c2*  (cos  (lx;U2*x)tu  >sh  (neta2*x)  )  t  , .  . 
r.  i  n  (beta2*>d  tsinh  (beta2*x) )  ; 

,/ 

1  1  'XV)  1  (1)“1; 

t  !  ncal  (2)  '-x; 

y  1  I  -rl oca  1  1  *  tempi*.-  t  l  ltr  dl  ’  *tenp7*n|ihl  *uniin'  ‘iifiii  '  *  tmi»>7  *  .  .  . 

117  hi  *iinti*n'  ‘mil*  '  *t  etip7*r  l.x  sl  ; 


i  * 

\  Hi  1 3  function  creates  tiie  p>  ixJoet  of  matrices  tor  int eyiatlon  of  MJi  % 
i  coefficient  matrix.  ' 

%  % 


fm  11  t  J  ini  y22“aat  r lng22  (x,  cl ,  c2,  Petal ,  tot a2,  f  1 ,  f 2,  t  1,  f -1,  rj  ■hi ,  w,  .  . 
wi  1 -s  id,  r  local,  t  turpi ,  t«tp2,  tm,  t  >,  r  local,  w:  lot,  untkjt ,  ytswec,  t  tup/,  .  .  . 
t «  m>tl,  t  enp9,  tempi  0,  tifilit  2,  st  i  f ! 

iH'hi  I  ),  i )  —  f  1  *  (cl  *  (cos  (total  *x)  teosh  (tot  at  *x)  )  tsln  (tot  a  l  *  x)  t  .  .  . 


KvS 


si nh (betal*x) ) +f3* (e2* (cos (beta2*x) icosh (beta2*x) )  f . .  „ 
sin  <beta2*x) +sinh (beta2*x) ) ; 

irphi  (3, 2)»=f2*  (cl*  (cos  (betal*x)  +cosh  (betal*x) )  +sin  (betal*x)  >. . . 
slnli  (betai*x) )  +f<i*  (c2*  (cos  (beta2*x)  (cosh  (beta2*x) 
sin (beta2*x> +sinh (beta2*x) ) ; 

rlocal (1)=1; 
rlocal  (2)=x; 


y22-r local '  *tenp8*nphi+um'  *npni '  *tenp8*nrphi; 


%%%  1i  %%%%%%%%%%%%%%%%%%%%%%%%t %%%%%%%%%%%%%%%%%%%%«%%%%%%%%%%%%%%%%%%%%%%%%%%% 


%  % 

%  Iliis  function  creates  tlie  product  of  matrices  for  integration  of  FQ  % 

%  coefficient  matrix.  % 

%  % 


function  y.33-astring33  (x,cl,c2,tetal/beta2/  fl,  f2,  fl,  f4,nphi, w,  . . . 
um's id,  riocai,  tenpl,  tenp2, um,  tenp,  rlocal, wdot,  umdot,  gra1--  c,tenp7, . . . 
teg  >fi,  terp9,  taipl  0,  nphi  2,  st  i  f ) 

m ’hi  (3,  l)*»fl *  (cl*  (cos (betal*x)  ♦cosli (bet at *x) )  »sin U«tal*x)  i 
iinlt  (Ixit.al  *x) )  if  3*  (c2*  (cos  (lx?ta2*x)  tc  ish  (beta2*x) )  i .  . . 
sin  (Ix;ta2*x)  +sinh (beLa2*x)  ); 

ill'll  i  (3, 2)-f2*  (cl  *  (cos  (beta*  *x)  tc< .sh  (l.x?ta 1 *x) ) t sin (l m!  a  1  * x)  t .  . . 
n  Inh  (!  »*lal  *x)  1  *  fl*  (e2*  (c  t  :>  (I  *.;ta2*x)  tcosfi  ((*;(  a2*x) )  » .  . .  , 

sin  (I xH a2*x)  *.si  nh (ix-t  a2*x) )  ; 

ii. "  d  (1)-1.0; 
r  lc »  .i  1.  (2)  “x; 

y  M  1  1 1  oc  il  *  ‘  t  ein>9'*t  1  oca  1  - 1  loca  1  '  *  I  enp9*n|>lii  *i«n  -  2 *  r  1  < x  a  1  '  *  tenp)  0*  .  .  . 
i.f  hi  *i»ix Jot  ixn'  *ii|9il  '  *t  enp9* t  loc.il  -tvN  *i;|3il  ’  *(  aip9*n|  9ii  *i»n-2 *  .  .  . 

» un '  'mill'  ‘teuplO'nphl  iinkit » t  1  cx  .i  1  *  *t«ip'  ‘gravvec  igi  awix; '  *  t  nip*  .  .  . 
ill'll  i  ‘urn; 

u:  nmmmmmutmnunumummmmmmmumummnu 

%  % 

%  lids  function  creates  tin?  product  of  nut  rices  (or  lnt  egt  at.  Ion  of  ft)  % 

\  i  .  >'f f  icient;  nwtiix.  \ 

*  % 


fun  t  inti  y44-ast  i  i  ikj44  (x,  ol,c2.  ,ioet  .,2,  f  1 ,  f2,  f  1,  f  4  .  !>•  ^ii,  w,  . 
wtc  :  lit,  t  local,  tcnpl,  torf>2,  um,  loifi.  r  local ,  tail  :  ,  •(*'151? 
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taTp8,te5rp9f  templO,iTpiii2,stl  f) 

rrphl  (3, l)~fl*  (cl*  (cos  (total  *x)  tcosh  (total *x) )  isin  (total *x)  ■»•... 
r;i  nh  (betal*x) )  +  f3*  (c2*  (cos  (beta2*x)  +cosh  (beta2*x) )  t. .  „ 
sin  (beta2,'rx)  +sinh  (beta2*x) )  ; 

iiphi (3,2)*=f2* (cl* (cos (betal*x) +cosh (betal*x) ) +sin (betal *x) +. . . 
sinh (tetal*x) ) +  f4*  <c2* (cos (beta2*x) +cosh (tota2*x) )  +  . . . 
sin (beta2*x) +sinh (beta2*x) ) ; 

y<l  <l=n  phi  '*w*  *w*nphi; 


%  This  function  creates  the  product  of  matrices  for  integration  of  GM  % 

'b  coefficient  nvatrix.  » 


function  y55=astring55  (x,cl,c2,total,beta2,  fl,  f2,  f3,  M,nphi,w, . . . 
wiesid,  r  .local,  tempi,  tenp2,um, tenp,  rlocaljWtotjUiirlotjgravvecytenp?,  . . . 

I:  ?np0,  tenp9,  tenpl.O,  nphi2,  st  if ) 

"Till  (3, 1)  -fl*  (cl*  (cos  (betal *x)  tcosh  (betal  *x) )  +sin  <betal*x) +. . . 

;;inh (total *x) » +f3* (c2* (cos (beta2*x) +cosh (beta2*x) )  t. . . 
r,  i  n  (beta2*x)  +sinh  (beta2*x) )  ; 

ii t  !il  (3, 2)  -=f2*  (cl*  (cos  (total  *x)  tcosh  (betal *x) )  +s in  (betal *x)  t. . . 
n  i  nh  (total *x)  >  f f4«(c2*  (cos  <tota2*x)  tcosh  (beta2*x) )  t. . . 
s  i n (beta2*x)  +sinh  (beta2*x) ) ; 

yh5-2*iif)lii 1  *w'*wciot*nphi; 

1  % 

l  I  Ids  function  creates  the  product:  of  matrices  for  integration  of  KN  % 

%  coefficient  matrix.  „ 

1  l 


functlor,  y66»ast:riiKj66(x, cl, c2, total ftota2,  fl,f2,  f3,  M,nj>hl,w,  . . . 
wi  i  '*> id,  i  local,  t«ipl,t«ip2,  tut,  tenp,  rlocal,wctot,imrJnt ,grawec,  t«tp7, 
tug)!!,  tai|>9,  LenplG,nphl2,stlf) 


ii(  ’ii  (  i,  1 )  -f  1  ■*  (cl  *  (cos  (total  *h)  tcosh  (total  *x)  )  tsi  n  (tot  a  I  *y  )  f  .  . 
s  i 11,1  (i  *•(  tl  *x> )  t  f  l*  (c2 *  (cos  (tot  a2*x)  fcosti  (tot  a2*x) )  t .  . . 
sin  (tota2*x>  Is  S  i ih  (tota2*x) )  ; 


nphi.  (3, 2)  ~f2*  (cl*  (cos  (betal  *x)  3 cosh  <betal*x) )  is  In  (betal*x)  1 . . . 
si,,),  (betal*x)  1  if 4*  (c2*  (cos  <bet:a2.*xl  fcosh  (beta2*x) )  + . . . 
sin  (beta2*x)  tsinh  (beta2*x) ) ; 


y  66-1(1*1 '  *w'  *wresid*mphi; 


%  % 

%  This  function  creates  the  product  of  matrices  for  integration  of  KN  % 

%  coefficient  matrix.  % 

%  % 


%%%%%%%%!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


function  y77--astrlng77  (x, cl, c2, betal, beta2,  f],  f2,  f3,  f4,rrphi,w, . . . 
wrosid,  rlocal,tenpl,teup2,um,  temp,  rlocal,  wdot,umdot,grawec,  tenp7, . . . 
teipO,  tenp9,  templO,  nphi2,  sti  f ) 

rrphi  2 (3, l)=fl* <betalA2) * (cl* (-cos (betal*x) lcosh(betal*x) ) . . . 

-sin (betal*x) +sinh (betal*x) ) if3* (beta2A2) * (c2* . . . 

(-cos  (beta2*x)  (cosh  (beta2*x) )  -sin  (beta2*x)  +sinh (beta2*x) ) ; 

mj:d i i  2  (3, 2)==f2*  (betalA2)  *  (cl*  (-cos  (betal *x)  1 cosh (betal *x) )  . . . 

-sin  (betal*x)  +s.inh (betal*x) )  +f4*  (beta2A2)  *  (c2*  . . . 

(-cos (beta2*x) +cosh (beta2*x) ) -sin (beta2*x) +sinh (beta2*x) ) ; 


y77-iiphi2 '  *stif  *nphi2; 


%r.n%%%%%%%%%%%%w%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%mm%%%%%%%%%%%%%%u%* 


%  j  % 

%  'Hi  is  function  creates  the  product  of  matrices  for  integration  of  FN  % 
%  coefficient  matrix.  % 
%  * 


function  yB8“33tringBB  (x,  cl,  c2,betal,bet  a2,  f  1,  f2,  f  3,  f  4,ri(>hl ,  w,  . .  . 
wrosid,  r  local, tenpl,tenp2,un,tenp,  r  1  oca  1 , wdot , umdnt: , g r awK, t: utp7 , . . . 
tnup8,tefip9,  tempi  0,11(4)1 2,  st  it ) 

rtf  '1 1 1  (3, 1  i~f  1*  (cl*  (cos  (bet  al *x)  icosh (betal*x) )  tain  (betal *x)  * .  . . 
s i  nli  (betal *x) )  if  3*  (c2*  (cos  (beta2*x)  icosh (bet a2*x> }  1. . . 
sin  ((»t.a2*x) tsinh (tjeta2*x) ) ; 

irphl  (3,  2)  2*  (cl  *  (cos  (betal*x)  icosh  (betal *xl )  tsl  n  (Ixitul  *x)  1 .  . . 

s  i  nil  (tel  a  I  *x) !  1  f4H  (c2*  (cos  (l  et  a2*x)  Icosh  (beta2*x)  )  t .  .  . 
sin  (beta2*x>  tsinh  (l;eta2*x) ) 


r 1 1  k  a 1 ( 1 ) -1 ; 


riocal  (2)  *=x; 

1 

yC8=nphi '  *vj'  *gcawec-fiphi '  *w ' 


*wresid* riocal 


00 
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